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Abstract 



Accretion disks are ubiquitous in the universe. Although difficult to observe directly, their presence 
is often inferred from the unique signature they imprint on the spectra of the systems in which 
they are observed. In addition, many properties of accretion-disk systems that would be otherwise 
mysterious are easily accounted for by the presence of matter accreting (accumulating) onto a 
central object. Since the angular momentum of the infalling material is conserved, a disk naturally 
forms as a repository of angular momentum. Dissipation removes energy and angular momentum 
from the system and allows the disk to accrete. It is the energy lost in this process and ultimately 
converted to radiation that we observe. 

Understanding the mechanism that drives accretion has been the primary challenge in accre- 
tion disk theory. Turbulence provides a natural means of dissipation and the removal of angular 
momentum, but firmly establishing its presence in disks proved for many years to be difficult. The 
realization in the 1990s that a weak magnetic field will destabilize a disk and result in a vigorous 
turbulent transport of angular momentum has revolutionized the field. Much of accretion disk 
research now focuses on understanding the implications of this mechanism for astrophysical ob- 
servations. At the same time, the success of this mechanism depends upon a sufficient ionization 
level in the disk for the flow to be well-coupled to the magnetic field. Many disks, such as disks 
around young stars and disks in binary systems that are in quiescence, are too cold to be sufficiently 
ionized, and so efforts to establish the presence of turbulence in these disks continues. 

This dissertation focuses on several possible mechanisms for the turbulent transport of angular 
momentum in weakly-ionized accretion disks: gravitational instability, radial convection and vor- 
tices driving compressive motions. It appears that none of these mechanisms are very robust in 
driving accretion. A discussion is given, based on these results, as to the most promising directions 
to take in the search for a turbulent transport mechanism that does not require magnetic fields. 
Also discussed are the implications of assuming that no turbulent transport mechanism exists for 
weakly-ionized disks. 
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"Since we astronomers are priests of the highest God in regard to the book of nature, it benefits 
us to be thoughtful, not of the glory of our minds, but rather, above all else, of the glory of God.' 



- Johannes Kepler 
Soli Deo gloria 
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Introduction 



Accretion disks form around gravitating objects because the angular momentum of the infalling 
gas is conserved. In order for the gas to accrete, however, its angular momentum must be removed. 
Understanding the mechanism underlying this angular momentum transport is a long outstanding 
puzzle in accretion disk theory. Much progress has been made in recent decades through the 
realization that a weak magnetic field will destabilize the flow in an ionized disk and result in a 
turbulent transport of angular momentum. One of the key features of this mechanism, however, 
is that the gas in the disk must be sufficiently ionized to couple to the magnetic field. It is likely 
that there are portions of disks, and perhaps entire classes of disks, in which the ionization is too 
low for the gas to destabilize. The search for a turbulent transport mechanism in weakly-ionized 
disks continues, therefore, to this day, in an effort to place our understanding of the evolution of 
these disks on as firm a theoretical footing as that for ionized disks. My research, in collaboration 
with Charles Gammie, has focused on several possible mechanisms: self-gravity, radial convection 
and vortices driving compressive disturbances. This dissertation summarizes our main results and 
discusses their relevance to the question of what drives accretion in weakly- ionized disks. 

The purpose of this chapter is to describe all of these ideas in detail and put them in their 
astrophysical context. I begin in HI. II with an overview of accretion disks: the systems in which 
they are observed and their general properties. The importance of angular momentum transport 
for the evolution of disks is discussed in Hl,2[ followed by a more detailed discussion of angular 
momentum transport in both ionized and weakly-ionized disks f§ H1.3l and I1.4JI . I give a brief 
overview in HI. 51 of the model that is used throughout the dissertation for analytic and numerical 

studies. HI .61 looks ahead to Chapter El in which I summarize and chart a course for future work. 1 

1 Portions of this chapter will be published in the proceedings of the Workshop on Chondrites and the Protoplan- 
etary Disk, November 8-11, 2004, Kauai, Hawaii. 
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1.1 Accretion Disks 



An accretion disk is a roughly cylindrical distribution of matter in orbit around a gravitating 
central object. It is supported against the gravitational pull of its host object primarily by the 
centrifugal forces arising from the angular momentum of the orbiting material. This support is 
slightly compromised, however, by the presence of dissipation or the application of external torques. 
As a result, angular momentum is redistributed through the disk and some of the disk material 
falls onto the central object, i.e., it accretes. The gravitational potential energy lost during this 
process is typically converted into radiation, which is the basis for our observations of accretion 
disk systems. Although disks are rarely observed directly (i.e., by being resolved in a telescope 
image), they imprint a unique signature on the spectra of their host systems. In addition, the 
accretion-disk paradigm easily accounts for many properties of astrophysical systems that would 
be difficult to explain otherwise. 

The material in an ac cretion disk covers a wide range of density and temperature scales 



(Bal bus and Hawlev 



1998). In most cases, collisional mean free paths are extremely short com- 



pared to the length scales of interest, and mean times between collisions are short compared to 
the time scales of interest. Disks are therefore usually modeled as a continuous fluid, using the 
macroscopic equations of gas dynamics rather than the microscopic equations of kinetic theory. If 
the fluid is ionized, it is referred to as a plasma. 

Accretion disks are found in a variety of astrophysical settings, including compact binary sys- 
tems (with a white dwarf, black hole, or neutron star), active galactic nuclei (AGN), and young 
stars. AGN consist of a supermassive black hole (M ~ 1O 8 M ) surrounded by an accretion flow. 
The presence of a disk in AGN is inferred primarily from the large luminosities of these systems, 
luminosities that cannot be accounted for by stellar nuclear burning but are easily provided by the 
large gravitational energy of the compact object. A spectacular exception to this indirect verifi- 
cation is the system NGC4258, an AGN in which the disk is directly observed via maser emission 



(Watson and Wallin 



19941 ) . Low Mass and High Mass X-Ray Binary (LMXRB and HMXRB) 



systems consist of a neutron star (NS) or black hole (BH) accreting matter from its companion; 
Cataclysmic Variables (CV) are binary systems in which the accreting object is a white dwarf (see 
Figure II. 1J) . The variability observed in these systems can be accounted for by models in which 
instabilities in a disk surrounding the compact object give rise to episodic accretion. Young Stellar 
Objects (YSO) are pre-main-sequence stars with a circumstellar disk, also known as protoplanetary 



2 



Table 1.1. Example Accretion-Disk Systems 



System 


Type 


M 

Mq 


M 
M yr- 


i 


Lace 
Lq 


R c (cm) 


T(K) 


H 

Rc 




NGC 4258 


AGN 


4 x 10 7 


1 x 10~ 


-2 


1 x 10 10 


1 x 10 18 


7 x 10 2 


2 x 10- 


-3 


Sco X-l 


LMXRB-NS 


1 


1 x 10" 


-8 


3 x 10 4 


1 x 10 10 


3 x 10 5 


5 x 10" 


-2 


LMC X-3 


HMXRB-BH 


10 


1 x 10" 


-8 


1 x 10 4 


1 x 10 11 


7 x 10 4 


3 x 10' 


-2 


UX Uma 


CV 


0.5 


1 x 10- 


-8 


1 x 10 1 


1 x 10 10 


2 x 10 5 


6 x 10- 


-2 


TW Hydrae 


YSO 


0.7 


5 x 10- 


10 


1 x 10" 2 


1 x 10 13 


1 x 10 2 


2 x 10" 


-2 



disks since they are thought to be the sites of planet formatio n. The presence of a disk in these 



systems is confirmed in many cases by direct observation (e.g.. 



Burrows et al 



1996). 




Figure 1.1 Schematic of a cataclysmic variable system. 

Table 11.11 lists typical values for various physical properties of these systems. The mass of the 
accreting object is denoted by M in Table 1X771 in units of the solar mass (M© = 2.0 x 10 33 g), and 
the accretion rate by M. The accretion luminosity 



GMM 



(1.1) 



(where Ri n is the inner radius of the disk and where G = 6.673 x 10 8 cm 3 g 1 s 2 is the gravitation 
constant) is given in units of the solar luminosity (Lq = 3.9 x 10 33 ergs _1 ). R c is a characteristic 



3 



radius for the accretion disk, and T is a characteristic temperature. The final column in Table ITTI 
contains the ratio of the vertical scale height H to the local radius. Here 



H 



n 



K 



(1.2) 



where c s is the isothermal sound speed and 



GM 



(1.3) 



is the local Keplerian orbital frequency. The values for T and H in Table 11.11 are obtained from 
the standard a-disk model, a derivation of which is outlined in the following section, using a = 



0.01. Values for the other quantities in Table 11.1 



Mivoshi et al 



van Paradiie 



)3) 



1995, 



Gammie et al 



1996; UX Uma: 



1999; Sco X-l: 



Frank et al 



1981 



were obtained fro m the litera ture (NGC4258 : 



Vrtilek et al 



1991 



; LMCX-3: 



TW Hydrae: iMuzerolle et al 



Paczvhski 



2000 



1983, 



Wilner et al 



2003). The standard disk model assumes that the internal energy of the disk material is efficiently 



radiated from the surfaces of the disk. One implication of this assumption is that the disk is typically 
quite thin, as can be seen from the last column of Table ITTT1 In addition, if the mass of the disk is 
much less than the mass of the central star, the orbital frequency of the gas Q = Qk + O(Hfr) 2 , 
so thin, low-mass disks have a nearly-Keplerian rotation profile. 



1.2 Angular Momentum Transport and Disk Evolution 



The accretion process consists of a net inward transport of matter and a net outward transport of 
angular momentum; a small fraction of the matter carries angular momentum outward, enabling 
the bulk of the matter to accrete. This angular momentum transport can take place 1) internally 
via a local exchange of momentum between fluid elements at adjacent radii or 2) externally via 
a global mechanism such as the rem oval of angular momentum by a wind off the surface of the 



disk (e.g. 



Blandfor d and Pavne 



1982). The focus of this dissertation is on the former: the internal 



diffusion of angular momentum. While global mechanisms such as winds and jets are certain to 
play a role in many accretion systems, their operation likely depends in a complicated manner upon 
the details of each particular system. Standard disk modeling typically igno r es their effects and 



assumes that angular momentum is transported internally (e.g., 



1991 



Sterzik and Morfill 



1994; 



Narita et al 



1994; 



Stepinski 



1997 



Pringle 



1981 



Gammie 



Ruden and Pollack 



1999). Whether or not 



it is possible to isolate this aspect of disk evolution and get meaningful results is a question that 
can only be answered as more comprehensive disk models are developed. 

As an example of the importance of global effects, as well as some of the difficulties involved 
in modeling them, consider the torques from MHD winds (axisymmetric pressure-driven winds 
have zero torque). Outflows are widely observed from YSO, and it is likely that the outflows are 
magnetically driven. Outflows are more common in young, high-accretion-rate systems. Highly 
uncertain estimates for the mass loss rate suggest that about 10% of the accreted mass goes into 
the jet and associated outflow. What is even less certain is the amount of angular momentum in the 
outflows, and therefore the role that they play in the evolution of disks on la rge scales fas opposed 



to the disk at radii less than a few tenths of an AU) . Wind models exist (e.g. 



1982 



Shu et al.l l994. 



2000 



Blandfor d and Pavne 



Wardle and Konis 1 1993), but there are large gaps n our understanding. 



We do not know what the strength of the mean vertical magnetic field, which organizes the wind, 
ought to be, nor how that mean field is transported radiall y through the disk, nor how the wind 



K5nigl and Pudritzl d.2000). 



evolves in time. A nice summary of this situation is given in 

Molecular shear viscosity is a natural mechanism for coupling fluid elements locally and trans- 
porting angular momentum internally, but the large Reynolds numbers of astrophysical flows (due 
to the large length scales involved) imply molecular shear viscosities that are much too tiny to 
account for the observations. The coupling due to molecular viscosity is simply too weak to explain 
the rapid variability and accretion rates that are observed. For example, the outburst duration 
in CV ranges from 2 — 20 da ys, with the interval between outbursts ranging from tens of days to 



tens of years ((Warner 



1995). The timescale for viscous diffusion over a distance / due to molec- 



ular viscosity is l 2 /v m , where v m is the molecular (or kinema tic) shear vis cosity. Using a value 



10 5 cm s 1 and a characteristic distance / = 10 10 cm (see 



Balbus 



2003 and Table H3J yields 



a viscous timescale of about 10 7 years, which is orders of magnitude too large to account for the 
timescales of CV outbursts. 

Standard disk modeling circumvents this problem by assuming the presence of an enhanced 
"anomalous viscosity" due to turbulence. The large Reynolds numbers are used to advantage, since 
our experience with laboratory flows indicates that the onset of turbulence typically occurs above 
a critical Reynolds number. The assumption of turbulent flow, along with the picture of turbulent 
eddies exchanging momentum with one another to drive accretion, underlies the majority of the 
phenomenological disk modeling that is currently used to explain observations. 
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The construction of a standard model for disk evolution proceeds as follows. We will use a 
cylindrical coordinate system centered on the accreting object, with r, <\> and z the radial, azimuthal 
and vertical coordinates, respectively. We start with mass conservation: 

% = -\ d r{rpv r )-^pv z ) (1.4) 

where p is the mass density and we assume axisymmetry (d/dcp = 0), on average. Integrating this 
equation vertically through the disk gives 

<9£ 1 

— = drir'EVr) + text (1.5) 

at r 

where £ = J dzp is the surface density, v r is a vertical average of the radial velocity, and T, ex t is the 
difference of —pv z evaluated at the upper and lower surface of the disk. It includes infall onto the 
disk, mass loss in winds, and mass loss through photoevaporation. It is positive when mass flows 
into the disk, and negative when mass flows out. 

The problem now is to find the radial velocity v r , which we can do using angular momentum 
conservation: 

Here pi is evidently the local density of angular momentum, and the right hand side of the equation 
is the divergence of an angular momentum flux density. H r ^ is a component of the stress tensor, 
sometimes referred to as the shear stress, with dimensions of pressure; it is the flux density of <\> 
momentum in the r direction. Likewise IL^ is the flux density of (ft momentum in the z direction. 
Again integrating vertically, 

= -\^-y W r<$> + rEM) + r + lt ext . (1.7) 

Here 



W r4> 



= j dzU. rt p — rT,v r l (1.8) 



is the integrated shear stress, but with one piece of it, proportional to the radial mass flux, peeled 
off. In models which assume that angular momentum transport is due to turbulence, W r( p is referred 
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to as the "turbulent shear stress." The external torque 

r = -rn^K " ffi«t (1.9) 

is the angular momentum flux into the upper and lower surface of the disk with one piece, propor- 
tional to the mass flux into the disk, peeled off. r includes the effects of, e.g., MHD winds; it is 
positive when angular momentum flows into the disk and negative when angular momentum flows 
out. In a steady state (d/dt = 0), the condition W r( p > must be met for an outward transport of 
angular momentum and inward accretion. 

The mass and angular momentum conservation equations ()1.5|) and Q1.7[) can be combined into 
a single equation governing the evolution of thin, Keplerian disks (multiply the continuity equation 
by rv,p, subtract the angular momentum equation, solve for v r , substitute back into the continuity 
equation) : 

^ 2 1 ' Vw*) - £) + (i.io) 



dt r dr \rVLdr " Q / 
If we assume that the shear stress W r $ is due to an anomalous viscosity u, and that the external 
torques and mass loss/infall are negligible, equation Ijl.lOj) becomes the basic equation for standard 
disk modeling: 

dZ 3 d ( 1 a . , n , 



dt r dr \ rQ dr 

In a steady state one can show that the accretion rate (inward mass flux = — 27rT,rv r ) is given by 

M = 3ttEz/. (1.12) 

In addition to setting r and T, e xt to zero, standard disk theory usually sets v = ac s H, which 
parameterizes our ignorance of W r $. If one reasonably assumes that the turbulent stress (an off- 
diagonal component of the stress tensor) must be associated with a pressure (an isotropic, diagonal 
component of the stress tensor), then a < 1. Most disk evolution models take a = const., or 
allow it to assume a few discrete values. For a disk around a solar- type star with a temperature 
of 300 K at 1 AU, this yields M = 9.9 x 10~ 9 aS M yr" 1 , or ~ 10~ 8 vr" 1 for a = 0.0 1 and 



S = 10 2 gcm 2 , roughly consistent with observed accretion rates (e.g. iGullbring et al.l . Il998) 
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1.3 Angular Momentum Transport in Ionized Disks 

While phenomenological disk modeling can proceed along its merry way without a clear demon- 
stration of the onset of turbulence in accretion disk flows, recent progress towards a first-principles 
understanding of disk turbulence has raised the possibility that more physically-motivated disk 
models can be developed. The primary breakthrough in our understanding came with the real- 
ization that the presence of a weak magnetic field destabilizes a disk on a dynamical time scale, 



resulting in the onse t of magnetohydrodynamic 



angular momentum ( Balbus and Hawlev 



1991 



(MHD) turbu 



199 



Balbu; 



ence and a vigorous outward flux of 



2003). 



This section gives a summary of the essential physics of this instability, generally termed the 
magneto-rotational instability, or MRI. The importance of ionization for the successful operation of 
the MRI in driving turbulence is also discussed, which leads in the following section to an overview 
of the main question addressed by this dissertation: what drives accretion in disks that are too 
weakly ionized to be unstable to the MRI? 

1.3.1 Magneto- Rotational Instability 

The MRI grows directly through exchange of angular momentum between radia lly-separated fluid 



eleme nts. This can be understood with a simple mechanical analogy introduced by 



Balbus and Hawlev 



(1992) and illustrated in Figure O 

Imagine that two small masses orbit with frequency 0(r) about a third, massive body. The 
masses are coupled by a spring; the natural frequency of the spring-mass system is 7. If 7 3> O, 
the bodies behave like a perturbed harmonic oscillator. But if 7 is lowered until 7 ~ 0, the orbital 
motion of the bodies begins to influence the dynamics, and something interesting happens. The 
outer mass has higher angular momentum but lower orbital frequency. It is pulled forward in its 
orbit by the spring; its angular momentum increases, moves to a higher orbit, and lowers its orbital 
frequency further. This stretches the spring, increases the rate at which the outer body gains 
angular momentum, and a runaway ensues. The outer body heads outward, acquiring angular 
momentum from the inner body. The inner body moves in a mirror image of the outer body as it 
loses angular momentum and falls inward. 

There is an exact correspondence between the modes of the spring-mass model and the MRI. 
One can think of the masses as fluid elements and the spring as magnetic field. The field has 



characteristic frequency 7 ~ va/A, where A is the separation of the masses and va = B/yJ Airp is 
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Figure 1.2 Mechanical analogy for the MRI. Two masses, orbiting in the direction indicated by the 
arrow about a massive body (bottom of the frame), are connected by a spring. The outer mass has 
higher angular momentum and lower orbital frequency. The lower mass is pulled back in its orbit 
by the spring, reducing its angular momentum. It sinks to a lower orbit, where it orbits faster, 
stretching the string and increasing the torque. A runaway ensues. The masses may be thought of 
as fluid elements connected by a magnetic field. 

the Alfven speed. va/\ < f2 implies instability. 

The simplicity of the dynamics shown in Figure 1131 suggests that the MRI is robust. Instability 
requires the presence of a weak (subthermal, i.e. B 2 /(8Tip) < c 2 ) magnetic field. A stronger 
magnetic field seems unlikely (it would likely be ejected from the disk by magnetic buoyancy), but 
if it were present it would likely be associated with other, even more powerful instabilities. The 
MRI also requires an angular velocity that decreases outward (dfl 2 /dlnr < 0), which is always 
satisfied in Keplerian disks (f2 oc r - 3 / 2 ). The maximum growth rate is ~ f2, independent of the 
field strength (unless diffusive effects are present). This surprising fact is easily understood once 
one realizes that the scale A of the instability decreases with the field strength: A ~ va/Q- 

What does the MRI tell us about W r( p, the key quantity in the evolution equation? Numeri- 
cal integration of the compressible MHD equations shows that the linear MRI initiates nonlinear 
turbulence. In the turbulent state one can measure the average value of 



where Sv^ is the noncircular azimuthal component of the velocity. There are two distinct contribu- 




(1.13) 
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tors to the shear stress: hydro dynamic velocity fluctuations, and magnetic field fluctuations, some- 
times referred to as the Reynolds and Maxwell stresses. Then a = const, x W r $j (£cf ) (the constant 



is a matt er of convention and takes several values in the literature). Loca" 



2000 



a 0.01 ( Hawlev et al 



1995 



Matsumoto and Taiima 



1995 



Ha wlev et al 



1997) . Global numer i cal models yield similar but slig htly larger values (e.g 



Machida et al 



2000 



Arlt and Riidiger 



2001). 



numerical models yield 



1996 



Matsuzaki et al 



Armit 



itaJ 199 fit 



Hawlev 



The presence of turbulence does not guarantee the W r< ^ > necessary for outward angular mo- 
mentum transport and accretion. For example, the turbu lence associated with vertical convection 



can produce a < (jStone and Balbu 



1996; 



Cabot 



1996]). Also, vortices in nearly incompressible 



disks produce a « (see Chapter EJ. The fact that MRI-driven turbulence has W r< j, > is thus 
a nontrivial result, although one that might have been anticipated because of the central role of 
angular momentum exchange in driving the linear MRI. 



1.3.2 MRI in Low-Ionization Disks 

MRI-generated MHD turbulence is the likely angular momentum transport mechanism in AGN and 



in binary systems during outbur sts. In portions of YSO disks, however, as wel 



X-Ray transients in quiescence (|Stone et al 



2000 



Gammie and Me nou , 



1998; 



as in CV disks and 



Menou . 



2000), the 



plasma is cool and nearly neutral. The conductivity of the gas is small by astrophysical standards 
and the field is no longer frozen into the gas. In some regions the field may be completely decoupled 
from the fluid, just as the Earth's lower atmosphere is decoupled from the Earth's magnetic field. 
Low ionization levels change the field evolution through three separate eff ects: Ohmic diffusion 



Hall d rift, and ambipolar diff u sion. Follow ing the beautiful discussion of 



200l|) (see also 



Wardle . 



19991 : 



Dcsch 



Balb us and Terquem 



2004), a measure of the correction to the field evolution 



equation (jl.21(l due to Ohmic diffusion is given by the magnetic Reynolds number RejJ = r//(vAH), 
where 77 = c 2 / (47r<7 e ) is the resistivity, c is the speed of light and the conductivity a e is proportional 
to the collision timescale for electrons with neutrals. Then 



Re M ^ 2 x 10 19 B 



r \ 3 / 2 (n 
AuJ 



M 



-1/2 



n 



-1/2 



(1.14) 



Here n is the neutral number density and n e is the electron number density. Ohmic diffusion 
destroys flux (via reconnection) and converts magnetic energy to thermal energy. 

Hall drift can be thought of as arising from the relative mean motion of the electrons and ions. 
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The associated correction to the field evolution equation is, for conditions appropriate to circum- 
stellar disks, typically comparable to Ohmic diffusion. Hall drift does not change the magnetic 
energy. 

Ambipolar diffusion arises from the relative mean motion of the ions and neutrals. Unlike Hall 
drift, it converts magnetic energy to thermal energy. The ratio of the ambipolar to Ohmic term 
~ 5 x 10 28 i? 2 T~ 1 / 2 n~ 2 (all cgs unit s; we have assumed that the number densities of electrons and 



ions are equal, but see 



Dcsch 



2004). In low-density environments such as the clouds from which 



young stars condense, ambipolar diffusion is the dominant nonideal effect; one can then think of 
the field as being locked into the ion-electron fluid, which gradually diffuses through the neutrals. 
At higher densities ambipolar diffusion becomes less dominant, and at the highest densities found 
in disks Ohmic diffusion dominates. Evidently the precise variation of the relative importance of 
these effects with location depends on the variation of ionization fraction, temperature, and field 
strength. In YSO disks, ambipolar diffusion tends to dominate in the disk atmosphere and at 
r > 10 AU. 

The linear theory of the MRI with Ohmic diffusion was first considered bv l.Tinl (|1996l ). Stability 



is recovered when the Ohmic diffusion rate rj/X 2 exceeds the growth rate of the instability va/X. 
Since A < H, this occurs when t]/{vaH) ~ 1. One can avoid expressing the stability condition in 
terms of the unknown field strength B, which likely arises through dynamo action induced by the 
MRI, by noting that for a subthermal field c s > va- Then when rj/(c s H) > 1 the disk is stable, 



lar disks 



although the MRI may be suppressed at even lower r\. Circumste 
are therefore stable to the MRI, over a large range in radius (e.g. 

The linear theory of the MRI with a mbipolar diffusion was first treated by 



Gammic 



rave j ]l (c K H) > 1, and 



1996). 



2003); 



1994) , and recon s idered more recently bv lDesch (2004); Ku nz and Balbud (|2004); 



Blaes and Balbus 



Salmeron and Wardlc 



Salmeronl (J200J). The natural expectation is that the MRI develops even in the presence of 



ambipolar diffusion as long as a neutral particle manages to collide with an ion (which can see the 
magnetic field ) at least once per orbit. A ssuming common values for the collision strengths and 



ion mass (e.g. 



Balb us and Terauem 



200 lh . this condition becomes 



A = 0.01 n e (r/AUf^iM/Mo)- 1 / 2 > 1. 



(1.15) 



The more recent round of papers points out that even when A < 1 there are unst able perturbations 



within a band of wavevectors k outside the purely axial wavevectors considered by 



Blaes and Balbus 



11 



(199J). These perturbations evade ambip qlar damping by orienting their magnetic field perturba- 



tions perpendicular to both B and k (see 



Desch 



2004). Instability can thus survive when A < 1, 



albeit in a narrow band of wavevectors and at greatly reduced growt h rates. 



The linear t heory of the MRI with Hall drift has been considered by 



Wardle! (| 19991 ) 



Balbus and Terauem 



(20011): iDeschl ((2004). There are always perturbations that become more unstable as Hall drift is 



turned on. The maximum growth rate is not affected. The combination of all these nonideal effects, 



togeth e r with a best gu e ss for the dis k ionization structure is discussed in 



(20o3); Sahncron (200 1): 



Dcsch 



(2004) 



Early numerical experiments by 



Safmeron and Wa rdle 



Haw lev et al 



( 19951 ) using a scalar resistivity, no explicit 



viscosity, no Hall drift and no ambipolar diffusion s uggested that t he M RI dynamo fails when 



c. q H/n < 10 4 . A simila r but more thorough study by 



Fleming et al 



(2000) found similar results. 



Sano and Stone! (2002a b) considered models that incorporated Ohmic diffusion and Hall drift, but 
not ambipolar diffusion (relevant in some regions of the disk). The most relevant of Sano &: Stone's 
mod els are probably t hose w ith net toroidal field or zero net field. Their results (see Figs. 14 and 



19 of 



Sano and Stone . 



2002b) suggest that Ohmic diffusion is the governing nonideal effect; a drops 
sharply when c s H/r] < 3 x 10 3 , and is only weakly dependent on the Hall parameter. 



1.4 Angular Momentum Transport in Weakly-Ionized Disks 

Accretion rates inferred from observations of weakly-ionized disks indicate that an enhanced vis- 
cosity or some other mechanism for angular momentum transport is operating in these disks. As 
discussed in the previous section, however, these same disks are likely to be stable to the MRI. This 
leaves open the question of what generates turbulence in weakly-ionized disks. As long as there was 
no firm theoretical understanding of the onset of turbulence in disks (ionized or not), it was reason- 
able to assume that all disks are turbulent due to their large Reynolds numbers. Laboratory shear 
flows are turbulent above a critical Reynolds number even though they are stable to infinitesimal 
perturbations (i.e., there is no linear instability to trigger the turbulence), and the extrapolation to 
astrophysical shear flows was a natural one to make. With the establishment of a robust transport 
mechanism in MRI-induced MHD turbulence, this assumption has come under critical scrutiny. If 
a mechanism can be established from first principles for ionized disks, it seems reasonable to main- 
tain the same standard for disks which are too weakly ionized to be MHD-turbulent. Although 
many attempts have been made, to date no robust transport mechanism akin to the MRI has 
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been established for low-ionization disks. This dissertation investigates three possible mechanisms 
in detail: gravitational instability, convection, and vortices driving compressive motions. Each of 
these mechanisms is summarized briefly here and discussed in detail in the remainder of the disser- 
tation. Since there are those who continue to argue for turbulence in disks by way of analogy with 
laboratory shear flows, a brief overview is also given of the current state of this controversy. 

1.4.1 Gravitational Instability 

Gravitational instability arises when self-gravity in the disk overcomes the stabilizing influences 
of pressure and rotation. The nonlinear outcome of this instability is either a gravito-turbulent 
state of marginal stability or fragmentation of the fluid into bound clumps. If a sustained gravito- 
turbulent state can be established, then steady outward angular momentum transport ensues. 
Instability tends to form temporary spiral enhancements in the density with a trailing orientation, 
and gravity then carries angular momentum along the spiral (there is a new term in W r s, a "New- 
ton" stress given by J dz g r g < f > /(8irG), where g r and are compon ents of th e grav itational field). 



Local numerical experiments exhibit a gravitational a up to ~ 1 (jGammie . 



20Mh . If this state 



can be maintained steadily through out the disk, it w ould provide an effective turbulent transport 



mechanism in weakly-ionized disks ()Paczvhski . 



1978). 



A sustained gravito-turbulent state cannot be established, however, if the cooling (due to radia- 
tion from the surfaces of the disk) is too strong. Then the clumps of matter formed by the instability 
cool before they can collide and heat each other via shocks. Fragmentation- the formation of small, 
bound clumps- results. The mean cooling time can be used to distinguish a gravito-turbulent disk 
from a fragmenting disk: 

r„ = jg. (1.16) 

where U = J dzu and u is the internal energy per unit volume, and A is the cooling function. The 
brackets () indicate an average over space and time, since the disk may have nonuniform density 
and temperature. 

Chapter |21 discusses in detail local numerical experiments which show that fragmentation occurs 
when r c fl < l. 2 These experiments also show that fragmentation occurs for a wide range of 
parameters, indicating that a gravito-turbulent state is difficult to sustain. In addition, cooling 
2 This result has bee n been demonstrated in other numerical experiments as well, both local and global (iGammieL 

EiiiilMHico ot „!lEnol;, 
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typically becomes more efficient with an increase in disk radius, making an extended, marginally- 
stable region unlikely. Gravitational instability thus does not appear to be a likely candidate for a 
turbulent transport mechanism in weakly-ionized disks. 



1.4.2 Convection 



As mentioned in M1.3.1[ vertical convection appears to transport angular momentum inward, op- 
posite to what is usually required for accretion. Indeed, arguments have been made that any 
incompressive disturbance or incompressible turbulen ce (of which convec tion is just one example) 



will drive angular momentum in the wrong direction (Balbu 



2000, 



2003). 



The possible role of radial c onvection in driving angular m ome ntum t r anspo rt has come to 



the fore recently with the work of 



Klahr and Bodenheimerl ((,2003) and 



Klahrl (|2004) on the "Global 



Baroclinic Instability" . One would expect that the combination of weak radial gradients and strong 
Keplerian shear in circumstella r disks would preclude any instabilities due to radial convection, yet 
Klahr and Bodenheimerl (|2003|) found turbulence and angular momentum transport in global hydrq - 



dynamic simulations with a modest radial equilibrium entropy gradient. The claim in lKlahrl ((2004) 
is that this activity, which grows on a dynamical time scale, is the result of a local hydrodynamic 
instability due to the presence of the global entropy gradient. 

Chapters|21and0]describe analytic and numerical work in a local model that attempts to confirm 
or refute these unexpected results. Chapter|3]describes a local stability analysis in radially-stratified 
disks, an analysis which uncovers no exponentially-growing instabilities for disks with a Keplerian 
rotation profile. Chapter |I] describes local numerical experiments which attempt to uncover any 
nonlinear instabilities that may be present. Disks with Keplerian shear are again found to be stable- 



It app ears, therefore, that the "Global Baroclinic Instability" claimed by 
([20031 ) is either global or nonexistent. 



Klahr and Bodenheimer 



1.4.3 Vortices 

The absence of a robust instability mechanism for generating hydrodynamic turbulence does not 
necessarily imply the absence of internal angular momentum transport. Chapter [5] describes nu- 
merical experiments which show significant shear stresses associated with finite-amplitude vortices 
that emit compressive waves and shocks. In these experiments, an initial field of random veloc- 
ity perturbations with Mach number ~ 1 forms anticyclonic vortices that provide an outward 
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flux of angular momentum corresponding to an initial a ~ 0.001 and decaying as t 1 / 2 . These 

results were obtained in a two-dimensional local model, and are likely to be modified consider- 

1 1 

ably by three-dimensional in stabilities, which tend to destroy two-dimensional vortices (Kerswell, 



2002 



Barranco and Marcu; 



2005). In addition, these results leave open the key question of what 



generates the initial vorticity. Both of these issues are discussed in more detail in Chapter 



1.4.4 Nonlinear Hydrodynamic Instability 

For some, the lack of a well-established mechanism for generating turbulence in weakly-io nized 



Richard and Z aim 



1999). In 



disks does not necessarily imply the absence of turbulence (e.g., 
the first place, our understanding of the onset of turbulence in simple laboratory shear flows is 
still incomplete, despite over a century of theoretical effort. Even when linear theory predicts 
stability of these flows at all values of the Reynolds number, experiments consistently show the 
onset of turbulence above a critical Reynolds number. The failure of linear theory to predict 
the outcome of experiments indicates that nonlinear instabilities (i.e., instabilities due to finite- 
amplitude disturbances) are the likely source of turbulence in these flows. Perhaps an analogous 
mechanism operates in weakly- ionized disks. In addition, since nonlinear stability is extremely 
difficult to prove due to the complexity of nonlinear dynamics, the question of stability in weakly- 
ionized disks remains, in some sense, an open question. As discussed in this section, however, 
no nonlinear instability mechanism has yet been established for a Keplerian shear flow, despite 
its apparent similarities with laboratory shear flows. In addition, the two main features that 
distinguish an accretion disk flow from a laboratory shear flow- rotational effects and the absence 
of rigid boundaries- seem to argue for the nonlinear stability of the former. 

The laboratory flow that most closely resembles an accretion disk flow is Couette- Taylor flow, 
which is flow between two concentric cylinders. Early t heoretical and experim ental results for this 
flow were obtained by Rayleigh, Couette and Taylor (see 



Dra zin and R .eid 



1981). Its linear stability 



is governed by the Rayleigh stability criterion, which states that a necessary and sufficient criterion 
for stability to axisymmetric disturbances is that 



^i^n(rf>o, 



(1.17) 



where k 2 is the square of the epicyclic frequency (also known as the Rayleigh discriminant). For a 
Rayleigh-stable flow, fluid elements displaced from circular orbits will undergo epicycles about their 
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equilibrium velocity at a frequency k. The stability criterion (|1.17f) is equivalent to the requirement 
that the specific angular momentum of the mean flow decrease with radius. 
Another la 



i boratory flow tha t has b een studied in depth is planar shear flow, also known as plane 
Couette flow ( Drazin and Reidl 



1981). This is flow between two parallel walls, the laminar state 



of which is a streamwise (parallel to the walls) velocity that varies linearly with distance from the 

walls. A necessary condition for the linear instability of a parallel shear flow with an arbitrary shear 

profile (i.e., an equilibrium velocity that is an arbitrary function of the coordinate perpendicular to 

the direction of flow) is that there be an inflection point in the equilibrium velocity profile. Since a 

linear shear profile does not meet this condition, plane Couette flow is also predicted to be stable 

based upon linear theory. 

Both Couette- Taylor flow and plane Couette flow show the onset of turbulence above a critical 

Reynolds number, against the predictions of linear theory. Theoretical efforts to explain this 

transition to turbulence have focused on 1) the transient amplification o f linear disturbances couple d 

I IF1 

with a nonlinear feedback mechanism to close the amplifier loop (e.g., Ba ggett and Trefethenl ll997); 
2) self-sustaining nonlinear proces ses that are t riggered at finite amplitude and are therefore not 



treatable by a linear analysis (e.g., Waleffe 
and secondary l inear instabi l ities ( e. g., iFarro ll ; nd loam ion 



1997) ; or 3) some comb ination of nonlinear mechanisms 



can be found in 



Bavlv et al 



(1988), 



Grossmann 



(2000) and 



1993); R eviews of these mechanisms 



Rempferl (|2003l ) . All of them include 



some aspects of the nonlinear dynamics and are generically referred to as nonlinear instabilities. 
While a discussion of their detailed operation is not necessary for the purposes of this dissertation, 
it is important to note that none of them has provided a complete understanding of the transition 
to turbulence in laboratory shear flows. 

T he app l icatio n of these ideas to accretion disks has continued sin c e the discovery of the 



MRI (Zahn, 



2001 



1991; 



Ri. 



Afshordi ct al 



20oi 



Dubrullc and Knobloch, 



2001a hT~garettI 



200 



4; 



Hersant et al 



2002 



1992, 



Mukhopadhvav et al 



200.4 



19q3: 



Dubrulle , 



Chagclishvili et al 



2004 



2003 



Richard and Davis 



Mu khopadhvav et al 



1200.4 



1993; 



Ioannou and Kakouri; 



200 



Umurhan et al 



; Richard, 2003; Klah 



r 



2004; 



2004; 



Umurhan and Regev 



2005). Much of this work 



has focused on the mechanism of transient amplification of linear disturbances coupled with non- 
linear feedback, since there are local nonaxisymmetric vortical perturbations which can experience 
an arbitrary amount of tran sient growth at infinite Reynolds number (a result that was recognized 



as early as 1907 by Orr; see 



Shepher 



1985). These solutions are discussed in detail in Chapter 01 
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where it is shown ( A3.4.3JI that an isotropic superposition of these perturbations has an energy 
that is constant with time. This seems to indicate that any potential mechanism for the onset of 
hydrodynamic turbulence in disks would be an entirely nonlinear process. Only nonlinear simu- 
lations can fully answer this question, however, and a full investigation of the effects of transient 
amplification over a wide range of initial perturbation amplitudes and spectra has not yet been 
made. To date, however, no numerical simulations have demonstrated a transition to turbulence 
from infinitesimal perturbations, and the results of 83.4.31 indicate that such a transition may not 
occur for a physically-realistic set of low-amplitude perturbations; 



1995, 



1229) found 



Early simulations of MHD turbulence in MRI-unstable disks (Ha wlev et al. 
that 1) when the magnetic fields were turned off the turbulence decayed away and 2) when ro- 
tational effects were removed, thereby converting the Keplerian flow into plane Couette flow, the 
turbulence increased and the magnetic field decayed away. While advocates of nonlinear instabil- 
ities in disk flows will o ften attribute th e absence of turbulence in hydrodynamic simulations to 



numerical diffusion (e.g. 



Longarett j[2002j), this l atter re s ult co nfirms the ability of these simulations 



to identify a nonlinear instability. In addition, 



Balbuf 



(2004j) has argued that due to a nonlinear 



scale invariance of the equations governing the local disk flow, any local instabilities that are present 
should be present at all scales and therefore not require high resolutions for their manifestation in 
local numerical simulations. 

Two subsequent comprehensive studi es of nonlinear instabilities in local numerical simulations 



(Balbus et al 



iquent 
■ Il99fi 



Ha wlev et al 



1999) have confirmed these results. Both Keplerian and plane 



Couette flows were investigated, using codes with very different diffusive properties, and rotational 
effects were cited as the key stabilizing factor in disks. One of the mechanisms for nonlinear 
instability in plane Couette flow is the generation of streamwise vortices by the shear, resulting 
in a secondary instability due to inflections in the spanwise (across the mean flow) direction. 
The epicyclic motions of fluid elements in a rotating flow prevent these streamwise vortices from 
developing. When rotational effects are removed, the nonlinear instability of planar shear flow is 
readily recovered. 

Before the work of 



Ha wlev et al 



(19993), the only Couette- Taylor (rotating- flow) experiments 
that showed the onset of turbulence had shear profiles that were sufficiently non-Keplerian to more 
closely resemble plane Couette flow than a rotationally-dominated flow (the shear profiles were near 
to linear instability, expression |1.17j ). More recent experiments, however, have shown the onset of 
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rd. 



2001a 



b), thus indicating 



turbulence in Couette- Taylor flow with a Keplerian shear profile (Richarc 
the presence of a nonlinear instability. These results, however, may simply highlight another key 
difference between labora tory shear flows and disk fl ows, namely the presence or absence of rigid 



boundaries. As noted in 



Garaud and Qgilvid (|2005j), early Couette- Taylor experiments revealed 



the importance of end effects in disturbing the laminar flow. Torque measurements in a system 
with an aspect ratio of 23 (the ratio of the cylinder lengths to the width of the gap between the 
cylinders) and an end plate corotating with the outer cylinder were 100% larger than measurements 
with the end plate stationary. An aspect ratio >, 40 was required to minimize the end effects. The 
ex perimental setup desc r ibed in Richard! ( 2001b ) has an aspect ratio of 25. 



Garaud and Qgilviel (|2005h proposed a model for the nonlinear dynamics of turbulent shear 
flows and also used their model to predict the onset of linear and nonlinear instability in shear 
flows both with and without rotation. The model accounts for many aspects of laboratory shear 
flow experiments. For reasonable model parameters, the model predicts nonlinear stability for 
Keplerian shear flows in the absence of boundaries and nonlinear instability for a wall-bounded 

y large Reynolds numbers. This is another 



" s " ,, " i """" 1 

Richard 



indication that the results observed by 



(2001b) may be due to boundary effects. 



1.5 Local Model 

Since the focus of this dissertation is on local mechanisms for angular momentum transport, all the 
analytic and numerical results are obtained in a local model of an accretion disk. Such a model 
can be obtained by a rigorous expansion of the fluid equations in |£c|/r, where x = (x,y,z) = 
(r — R ,R ((j) — 4> — Q.(R )t),z) ~ 0(H) are the local Cartesian coordinates of the fluid with 
respect to a fiducial radius R Q and fiducial angle 4> Q (see Figure H~3|) . Since the local coordinates are 
assumed to vary on the order of the disk scale height H, the local model expansion is only valid for 
thin disks with H/r <C 1 (see Table ITT]) . This local frame is corotating with the fluid in the disk at 
a distance R Q from the central object and at a frequency Q(R ), the local rotation frequency of the 
disk. Local curvature is neglected, but centrifugal and Coriolis forces are retained. The additional 
simplifying assumption of an infinitesimally thin disk is made, which implies a vertical integration 
of the fluid variables. 



18 




Figure 1.3 Local coordinate system at t = 0. 

The resulting equations of motion for a fluid in the local model (including self-gravity) are 

1 / (B V).B 

d t v + (v ■ V)v = V IP + — + - — - — - - 2Slxv + 3n 2 xx, (1.18) 

S \ 8ir J Airp 

where v is the fluid velocity with respect to the rotating frame, B is the magnetic field, P is the 
two-dimensional pressure, X is the column density and 4> is the disk potential with the time-steady 
axisymmetric component removed. The last two terms on the right-hand side of equation (|1.18|) 
incorporate the effects of Coriolis and centrifugal forces as well as the gravitational acceleration 
due to the central point mass and the time-steady axisymmetric component of the disk. These 
equations of motions are valid for a disk system in which the gravitational potential is dominated 
by the central object; the fluid in such a disk follows a Keplerian rotation curve, iu ~ r -1 / 2 . 
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The continuity, internal energy and induction equations retain their usual form: 



<%E + 0- V)£ + £V-r> = 0, 



(1.19) 



d t U + (v ■ V)U +(U + P) V • v + A = 0, 
(where £/ is the internal energy per unit area and A is the cooling function) and 



(1.20) 



d f B - Vx(vxB) = 0. 



X21) 



Equations (|1.18j) through (|1.21jl are the equations of compressible ideal magnetohydrodynamics 
(MHD) in the local model. Magnetic fields are assumed to be dynamically unimportant for most 
of research described in this dissertation, in which case equations (|1.18|) through (|1.20|) reduce 
to the equations of compressible hydrodynamics. The disk self-gravity ((f) in equation |1.18| ) and 
explicit cooling (A in equation jl.20| ) are neglected except in the work discussed in Chapter |2j In 
Chapters HI - EJ the fluid is assumed to obey an ideal-gas equation of state, P = (7 — 1)U, where 7 
is the adiabatic index. Chapter \5\ assumes an isothermal equation of state P = c^S. 

With constant density and pressure in equilibrium, an exact steady-state solution to equations 
(jl,18|) through (|1.20j) is v Q = —^Qxy. This uniform shear velocity is a manifestation of differential 
rotation of the fluid in the disk. As a result, the (two-dimensional) local model is referred to as the 
"shearing sheet" . The numerical implementation of the shearing sheet requires a careful treatment 
of the b oundary condi t ions i n the radial direction. These boundary conditions are described in 



detail in 



Haw lev et al 



(1995j). In brief, one uses strictly periodic boundary conditions in y and 
shearing-periodic boundary conditions in x. The latter is done by enforcing periodic boundary 
conditions in the radial direction followed by an advection of the boundary fluid due to the shear. 
This assumes that the shearing sheet is surrounded by identical boxes that are strictly periodic 
initially, with a large-scale shear flow present across all the boxes. 



The shearing-sheet equations are evolved using a ZEUS-based scheme ((Stone and Norman . 

1992! ): a time-explicit, operator-split, finite-difference method on a staggered grid which uses an 



artificial visco sity to captur e shocks. An important modification of the standard shearing sheet, 



introduced by 



Massetl ((20,00), is the splitting of the overall shear velocity from the rest of the flow. 



This overcomes a practical limitation of the standard shearing sheet, which is the small Courant- 
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limited time step imposed by the large shear velocities at the edges of the sheet; for numerical 
stability of grid-based schemes the Courant condition requires time steps to be lower than the grid 
spacing divided by the maximum velocity on the grid. The larger the box, the more severe this 
limitation becomes. Separating out the shear removes this limitation and allows one to increase 
the size of the shearing sheet arbitrarily. This separation is done by replacing v y by v + 5v y in 
the fluid equations, and then evolving 5v y ; this can be done because there is no evolution of v 
directly {dtv = and Vv D = constant). Advection of other fluid variables by v Q is done by splitting 
the distance over which the fluid is sheared into an integral and fractional number of grid zones: 
the fluid variables are simply shifted an integral number of zones and then advected in the usual 
manner for the remaining fractional part (which does not require a higher effective velocity than 
any other part of the flow). 



1.6 Discussion 



While a rigorous proof of the stability of weakly-ionized disks may well be impossible, the results of 
this dissertation add to the already strong evidence against a turbulent angular momentum trans- 
port mechanism in weakly-ionized disks. Gravitational instability likely results in fragmentation, 
radial convection is suppressed by differential rotation and two-dimensional vortices, which provide 
a decaying flux of angular momentum, are likely to be unstable in three dimensions. Evidence 
against a local, nonlinear, purely hydrodynamic instability is mounting; 



Accretion may be driven glo bally by a magneto-centrifugal wind 



1989; 



Livio and Spruit, 



or tidally- induced spiral waves ((Larson . 

excited by planets embedded in the disk jGoodman and Haliko\ 



Papaloizou and Pringk 



There also exist global instabilities (e.g., 

amount of turbulence and angular momentum transport (e.g., 
are instabilities associated with the dust layer in YSO disks (e.g., 



Blandfor d and Pavne 



1982) 



1991), or locally via spira l wave s 



2001 



1984 



Sari and GolHreichl. |2Q0J) . 
IB that result in a small 



Hawlev 



1987 ). In addit ion, there 



Garaud and Lin 



2004) that will 



generate some amount of turbulence in those systems. While one or more of these mechanisms 
may play a role in transporting angular momentum in certain systems, their dependence upon 
global structure or other special features in order to operate makes their broad application to 



weakly-ionized disks dou btful. A 



in ionized surface layers (Gammie, 



ternatively, weakly-ionized disks may simply be inactive except 



1996). 



A detailed discussion of these possibilities is beyond the scope of this dissertation, but a brief 
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discussion of layered accretion is given in Chapter along with some proposals for future modeling 
based upon that idea. Chapter also summarizes the main results and implications of this work, 
and provides some direction as to where to go from here in the search for a turbulent angular 
momentum transport mechanism in weakly-ionized accretion disks. In addition, proposals are 
made for future investigations of the properties of turbulent stresses in ionized disks, with a view 
towards incorporating these properties in advanced, physically-motivated disk models. 
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Nonlinear Outcome of Gravitational 
Instability in Disks with Realistic 
Cooling 



2.1 Chapter Overview 

We consider the nonlinear outcome of gravitational instability in optically-thick disks with a re- 
alistic cooling function. We use a numerical model that is local, razor-thin, and unmagnetized. 
External illumination is ignored. Cooling is calculated from a one-zone model using analytic fits 
to low temperature Rosseland mean opacities. The model has two parameters: the initial surface 
density S Q and the rotation frequency f2. We survey the parameter space and find: (1) The disk 
fragments when {(t c ))Q ~ 1, where ((t c )) is an effective cooling time defined as the average internal 
energy of the model divided by the average cooling rate. This is consistent with earlier results 
that used a simplified cooling function. (2) The initial cooling time r co for a uniform disk with 
Toomre stability parameter Q = 1 can differ by orders of magnitude from ((t c )) in the nonlin- 
ear outcome. The difference is caused by sharp variations in the opacity with temperature. The 
condition t co Q ~ 1 therefore does not necessarily indicate where fragmentation will occur. (3) 
The largest difference between ((r c )) and r co is near the opacity gap, where dust is absent and 
hydrogen is largely molecular. (4) In the limit of strong illumination the disk is isothermal; we find 
that an isothermal version of our model fragments for Q < 1.4. Finally, we discuss some physical 
processes not included in our model, and find that most are likely to make disks more susceptible 
to fragmentation. We conclude that disks with ((t c }}Q < 1 do not exist. 1 



2.2 Introduction 



The outer regions of accretion disks in both active galactic nuclei (A GN) and young stellar 



Shlos mau et al .11199(1 : YSOs: 



objects 



(YSO) are close to gravitational instability (for a review see, for AGN: 

1 Published in ApJ Volume 597, Issue 1, pp. 131-141. Reproduction for this dissertation is authorized by the 
copyright holder. 
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Adams and Linl ll993). Gravitational instability can be of central importance in disk evolu tion. In 



some disks, it leads to the efficient redistribution of mass and angular momentum (e.g 



1984 



Laughlin and Rozvczka 



1996 



Larson 



Gammiell2n0ll ). In other disks, gravitational instability leads to 



fragm entation and the formation of bound objects. This may cause t he truncation of circumnuclear 



disks (Goodman 



2003), or the formation of planets (e.g. 



Boss 



1997, and references therein). 



We will restrict attention to disks whose potential is dominated by the central object, and whose 
rotation curve is therefore approximately Keplerian. Gravitational instability to axisymmetric 
perturbations sets in when the sound speed c s , the rotation frequency O, and the surface density 
S satisfy 



Q 



<Q. 



crit 



1 



(2.1) 



( Toomrc 



1964; 



Goldreich and Lvnden-Belll 



1965). Here Q cr it = 1 for a "razor-thin" (two- dimen- 
sional) fluid dis k model of the sort we will consid er below, and Q cr it = 0.676 for a finite-thickness 
isothermal disk ( Goldreich and Lvnden-BelJ . 19651 ) . 2 The instability condition (|2.1j) can be rewrit- 



ten, for a disk with scale height H ~ c s /f2, around a central object of mass M*, 



M disk > -M*, 



(2.2) 



where M^sk = vrr 2 S. For YSO disks H/r ~ 0.1 and thus a massive disk is required for instability. 
AGN disks are expected to be much thinner. The instability condition can be rewritten in a third, 
useful form if we assume that the disk is in a steady state and its evolution is controlled by internal 
("viscous") transport of angular mom entum. Then the accret i on ra te M = 37ra:c 2 £/ri, where a < 1 



is the usual dimensionless viscosity of 



• . 3qc^ 
M > 



Shakura and Sunvaevl (|1973h . and 



G 



implies gravitational instability (e.g. 



7.1 x 10" 4 a 



1 km s 1 



M Q yr- 1 



(2.3) 



Shlosman et al 



(1990)). Disks dominated by external torques 



(e.g. a ma gnetohydrodyna mic [MHD] wind) can have higher accretion rates (but not arbitrarily 



higher; see 



Goodman 



2003) while avoiding gravitational instability. 



For a young, solar-mass star accreting from a disk with a = 10 2 at 10 6 M Q yr , equation 

(|2.3j) implies that instability occurs where the temperature drops below 17 K. Disks may not be 

2 Fo r global mo dels with radial structure, nonaxisymmetric instabilities typically set in for slightly larger values of 
Q (see Boss 1998 and references therein). 
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this cold if the star is located in a warm molecular cloud where the ambient temperature is greater 
than 17 K, or if the disk is bathed in scattered infrared light from the c entral star (althou gh there 
is some evidence for such low temperatures in the solar nebula, e.g. lOwen et al.N l999). If the 
vertically-averaged value of a is small and int ernal dissipation is confined to surface layers, as in 
the layered accretion model of iGammid (J1996J), then instability can occur at higher temperatures, 
although equation (|2.2[) still requires that the disk be massive. 

AGN disk heating is typically dominated by illumination from a central source. The temperature 
then depends on the shape of the disk. If the disk is flat or shadowed, however, and transport is 



4258 ( 


Mivoshi 


et al.. 


Gammie et al. 


1999) 



1995) the accretion rate may be as large as 10 2 M@yr 1 (Las ota et al 



1996; 



1999). Equation ((2~3)) then implies that instability sets in where T < 10 4 (a/10~ 2 ) K. 



If the disk is illuminatio n-dominated th en Q fluctuates with the luminosity of the central source. 



In a previous paper (jGammk 



2001J), one of us investigated the effect of gravitational instability 
in cooling, gaseous disks in a local model. A simplified cooling function A was employed in these 
simulations, with a fixed cooling time r co : 



A U 

A = , 

T~co 



(2.4) 



where U = the internal energy per unit area. Disk fragmentation was observed for r co il < 3. The 
purpose of this paper is to investigate gravitational instability in a local model with more realistic 
cooling. 

Several recent numerical experiments have included cooling, as op posed to isothermal or adia- 



batic evolution, and we can ask whether these results are consistent with 



Gammiel fcOQlh 



Nelson et al 



(2000) studied a global two-dimensional (thin) SPH model in which the vertical density and tem- 
perature structure is calculated self-consistently and each particle radiates as a blackbody at the 
surface of the disk. The initial conditions at a radius corresponding to the minimum initial value 
of Q (~ 1.5) for these simulations were E G ~ 50gcm~ 2 ,il ~ 8 x 10~ 10 s _1 ; the initial cooling time 
un der these circumstanc es is r co ~ 250 Q , so fragmentation is not expected and is not observed. 

(2001) consider a global three dimensional (3D) Eulerian hydrodynamics model 



Durisen et al 



in which the volumetric cooling rate varies with height above the midplane so as to preserve an 
isentropic vertical structure. The cooling time is fixed at each radius. Thei r cooling t ime > 10O — 1 
at all radii, so fragmentation is not expected based on the criterion of IGammid ([20011 ) . The 
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si mulations show str ucture formation due to gravitational instabilities but not fragmentation. 



Rice et al 



(2003) consider a global 3D SPH model with a cooling time that is a fixed multiple 
of (r). They find that their disk fragments when r co ~ 3il _1 and M^isk = 0.1M*. For a more 
massive disk {Mdisk = 0.25M*), fragmentation occurred at somewhat higher cooling tim es (r m ~ 
10T2" 1 ). This is effectively a global generalization of the local model problem considered bv lGammie 



(200 1f). The fact t hat the results are so consistent suggests that the local, thin approximation used 



m 



Gammie (2001 ) and here give a reasonable approximation to a global outcome. 



Maver et al 



(2002) consider a global three dimensional SPH model of a circumstellar disk. 
Explicit cooling is not included, but the equation of state switches from isothermal to adiabatic 
when gravitational instability begins to set in. This is designed to account for the inefficient cooling 
of dense, optically thick regions. Fragmentation is observed. Realistic cooling can have a complex 
influence on disk evolution, and it is not clear that switching between isothermal and adiabatic 



behavior "brackets" the outcomes that might 
Other notable recent work, such as that by 



je obtained when full cooling is used. 



Boss! (j2002j) , includes strong radiative heating in the 



sense that the effective temperature of the external radiation field Tj rr is comparable to or larger 
than the di sk midplane tem perature T c . In the limit that Tj rr <C T c we recover the limit considered 



here and in 



Gammid (|200lh : in the limit that Tj rr 3> T c the disk is effectively isothermal. 



The plan of this paper is as follows. In ^12.31 we describe the model, with a detailed description 
of the cooling function given in M2.4I The results of numerical experiments are described in < j2.5l 
Conclusions are given in ^12.61 



2.3 Model 



Gammiel (|200ll ) in every respect except that 



The model we use here is identical to that used in 
we use a more complicated cooling function. To make the description more self-contained, we 
summarize the basic equations of the model here. The model is local, in the sense that it considers 
a region of size L where L/r Q <C 1 and r Q is a fiducial radius. We use a local Cartesian coordinate 
system x = r — r Q and y = {(f) — £lt)r D , where r, (ft are the usual cylindrical coordinates and f2 is the 
orbital frequency at r Q . The model is also thin in the sense that matter is confined entirely to the 
plane of the disk. 

Using the local approximation one can perform a formal expansion of the equations of motion in 
the small parameter L/r Q . The resulting equations of motion read, where v is the velocity, P is the 



2G 



(two-dimensional) pressure, and is the gravitational potential with the time-steady axisymmetric 
component removed: 

r>ii up 

(2.5) 



Dv 
~Dt 



VP n 

■— 2ft x v + 3tt 2 xx - Vd>. 



For constant pressure and surface density, v = —^Qxy is an equilibrium solution to the equations 
of motion. This linear shear flow is the manifestation of differential rotation in the local model. 
The equation of state is 

P = (7 - 1)U, (2.6) 



where P is the two-dimensional pressure and U the two-dimensional internal energy. The two- 
dimensional (2D) adiabatic index 7 can be mapped to a 3D adiabatic index r in the low-frequency 
(static) limit. For a non-self-gravitating disk 7 = (3r — l)/(r + 1) (e.g. 



Qstriker et al 



Goldreich et al 



1986 



1992). For a strongly self-gravitating disk, one can show that 7 = 3 — 2/T. We 



adopt r = 7/5 throughout, which yields 7 = 11/7. 
The internal energy equation is 



dU 
~dt 



+ V • (Uv) 



-PV ■ v - A, 



(2.7) 



where A = A(S, U, Q) is the cooling function, fully described below. Notice that there is no heating 
term; heating is due solely to shocks. Numerically, entropy is increased by artificial viscosity in 
shocks. 

The gravitational potential is determined by the razor-thin disk Poisson equation: 



V 2 = 4vrGS5(z) 



(2i 



For a single Fourier component of the surface density this has the solution 



2irG 



S„ik-x—\kz 
k e 



(2.9) 



A finite-thickness disk has weaker self-gravity, but this does n ot qualitatively change the dynamics 



1965) 



of the disk in linear theory (Gol dreich and Lvnden-Belll 

We integrate the governing equations using a self-gravitating hydrodynamics code based on 



ZEUS (|Stone and Norman 



1992). ZEUS is a time-explicit, operator-split, finite-difference method 



on a staggered mesh. It uses an artificial viscosity to capture shocks. Our implementation has 
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been tested on standard linear and nonlinear problems, such as so und waves and shock tubes. We 
use the "shearing box" boundary conditions, described in detail by 



Ha wlev et al 



( 1995|), and solve 



the Poisson eq uation using the Fourier transform method, modified for the shearing box boundary 



Gammid (|200l|) for further details on boundary conditions, numerical methods and 



conditions. See 
tests. 

The numerical model is always integrated in a region of size L x L at a numerical resolution of 
N x N. In linear theory the disk is most responsive at the critical wavelength 2c^/GT> . 3 We have 
checked the dependence of the outcome on L and found that as long as L > 2c1/GT, a the outcome 
does not depend on L. We have also checked the dependence of the outcome on ./V and found that 
the outcome is insensitive to N, at least for the models with N > 256 that we use. 



2.4 Cooling Function 



Our cooling function is determined from a one-zone model for the vertical structure of the disk. 
The disk cools at a rate per unit area 

A = 2o-T e 4 , (2.10) 

which defines the effective temperature T e . The cooling function depends on the heat content of 
the disk and how that content is transported from the disk interior to the surface: by radiation, 
convection, or perhaps some more exotic form of turbulent tr ansport such as MHD waves. Low 



temp e rature d i sks ar e expected to be convectively unstable (e.g. 



1980). 



Cassenl (|1993T ) has argued, however, that the radiative heat flux in an adiabatically-stratified 



Lin and Papaloizou 



disk is comparable to the heat dissipated by turbulence (in an a-disk model), suggesting that 
convection is incapable of radically altering the vertical structure of the disk. We will consider only 
radiative transport. 



If the disk is optically thick in the Rosse 



treated in the diffusion approximation, then ([Hubenv 



and mean sen se, so that radiative transport can be 



1990) 



8T[ 
3 r 



(2.11) 



where r is the Rosseland mean optical depth and T c is the central temperature. We will assume 



The wavelength corresponding to the minimum in the dispersion relation for axisymmetric waves. 
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that T c ~ T, where 

T= ^A^ (2.12) 

and 

^ = 7(7"1)|, (2.13) 

which follows from the equation of state and the assumption that the radiation pressure is small 
(we have verified that this is never seriously violated). Here ks is Boltzmann's constant, m p is the 
proton mass, and p is the mean mass per particle, which we have set to 2.4 in models with initial 
temperature below the boundary between the grain-evaporation opacity and molecular opacity and 
p = 0.6 in models with initial temperature above the boundary. 
The optical depth is 

/•oo 

t= / dzK( Pz ,T z ) Pz (2.14) 







where k is the Rosseland mean opacity, p z and T z are local density and temperature, and z is the 
height above the midplane. Following the usual one-zone approximation, 

POO 

/ dzK(p z ,T z )p z a HK(p,f)p (2.15) 
Jo 

where the overbar indicates a suitable average and H ~ c s (T)/Q is the disk scale height (we ignore 
the effects of self-gravity on the disk scale height, which is valid when locally Q > 1). Taking T ~ T 
and p ~ Jj/(2H) then gives a final, closed expression for A. 



Bell and LiiJl|l994h . 



We have adopted the analytic approximations to the opacities provided by 
These opacities are dominated by, in order of increasing temperature: grains with ice mantles, 
grains without ice mantles, molecules, H~ scattering, bound- free/free- free absorption and electron 
scattering. The molecular opacity regime is commonly called the opacity gap; it is too hot for dust, 
but too cold for H~ scattering to contribute much opacity. The opacity can be as much as 4 orders 
of magnitude smaller than the ~ 5 gem -2 typical of the dust-dominated opacity regime. It turns 
out that this feature plays a significant role in the evolution of gravitationally-unstable disks. 

To sum up, the cooling function is 

ACE,U,n) = — — . (2.16) 
3 r 
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Table 2.1. Scaling Exponent for Cooling Time as a Function of Surface Density 



Opacity Regime 


a 


b 


Exponent 


Ice grains 


U 


o 
Z 


in/7 

1U/ 1 


Evaporation of ice grains 





-7 


-26/7 


Metal grains 





1/2 


4/7 


Evaporation of metal grains 


1 


-24 


-89/7 


Molecules 


2/3 


3 


52/21 


H scattering 


1/3 


10 


131/21 


Bound-free and free-free 


1 


-5/2 


-3/7 


Electron scattering 








2/7 



For a power-law opacity of the form k = Kop a T b , this implies that 

A ~ ^-^a/2+bjji+a/2-b ^ ^ ^ 

From this it follows that the cooling time r c = U/A scales as 

T c ~ s5+3a/2-6 J j-3-tt/2+6 i ( 2 .18) 

If the disk evolves quasi-adiabatically (as it does if the cooling time is long compared to the 
dynamical time) then U ~ £ 7 and 

t c ~ s 5 " 37+(a/2)(3 ~ 7)+fe(7 ~ 1) . (2.19) 

Table l2~T1 gives a list of values for this scaling exponent for our nominal value of 7 = 11/7. Notice 
that, when ice grains or metal grains are evaporating, and in the bound- free/free- free opacity 
regime, cooling time decreases as surface density increases. 

Our cooling function is valid in the limit of large optical depth (r 3> 1). Since the disk becomes 
optically thin at some locations in the course of a typical run, we must modify this result so that 
the cooling rate does not diverge at small optical depth. A modification that produces the correct 
asymptotic behavior is 

This interpolates smoothly between the optically-thick and optically-thin regimes and is propor- 
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tional to the (Rosseland mean) optical depth in the optically-thin limit. While it would be more 
physically sensible to use a Planck mean opacity in the optically-thin limit, usually the optically- 
thin regions contain little mass so their cooling is not energetically significant. An exception is in 
the opacity gap, where even high density regions become optically thin. 

Our simulations begin with S and U constant. The velocity field is perturbed from the 
equilibrium solution to initiate the gravitational instability. The initial velocities are v x = Sv x , 
v y = — |f2x + 5v y , where 5v is a Gaussian random field of amplitude (Sv 2 ) /c 2 = 0.1. The power 
spectrum of perturbations is white noise (v 2 . ~ k°) in a band in wavenumber k cr u/4: < \k\ < 4:k cr u 
surrounding the minimum k cr a = l/(irQ 2 ) (with G = S D = £2 = 1) in the density-wave dispersion 
relation. We have checked in particular cases that for 10~ 3 < (5v 2 ) /c 2 < 10 the outcome is qualita- 
tively unchanged. This is expected because disk perturbations (unlike cosmological perturbations) 
grow exponentially and the initial conditions are soon forgotten. 

Excluding the initial velocity field, the initial conditions for a spatially- uniform disk consist of 
three parameters: T, ,U , and f2. We fix Q = 1, leaving two d egrees of freedom . In models with 



Gammid (|200l|), these degrees of 



simple, scale-free cooling functions such as that considered by 
freedom remain and can be scaled away by setting G = S G = Q = 1. That is, there is a two- 
dimensional continuum of disks (with varying values of T, Q and f2, but the same value of Q) that 
are described by a single numerical model. 

The opacity contains definite physical scales in density and temperature. The realistic cooling 
function considered here therefore removes our freedom to rescale the disk surface density and 
rotation frequency. That is, there is now a one-to-one correspondence between disks with fixed X 
and O and our numerical models. 

The choice of S Q and $7 as labels for the parameter space is not unique. Internally in the code 
we fix the initial volume density (in g cm -3 ) and the initial temperature (in Kelvins). These choices 
are difficult to interpret, however, since they are tied to quantities that change over the course of 
the simulation; f2 and the mean value of X do not. 

The cooling is integrated explicitly using a first-order scheme. The timestep is modified to 
satisfy the Courant condition and to be less than a fixed fraction of the shortest cooling time on 
the grid. We have varied this fraction and shown that the results are insensitive to it, provided 
that it is sufficiently small. 
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2.5 Nonlinear Outcome 



2.5.1 Standard Run 

Consider the evolution of a single "standard" run, with S c = 1.4 x 10 5 gcm~ 2 and $7 = 1.1 x 
10~'sec . This corresponds to T Q = 1200 and r co = 9.0 x 10 4 $7 _1 . The model size is L = 
320G£ o /fi 2 and numerical resolution 1024 2 . The model initially lies at the lower edge of the 
opacity gap. 

The evolution of the kinetic, gravitational and thermal energy per unit area ({E^), (E g ) and 
(E t h) respectively) normalized to G 2 Sq/0 2 , 4 are shown in Figure l2~Tl After the initial phase of 
gravitational instability the model settles into a statistically-steady, gravito-turbulent state. It does 
not fragment. Cooling is balanced by shock heating. Energy for driving the shocks is extracted 
from the shear flow, and the mean shear flow is enforced by the boundary conditions. 



t — i — i — i — | — i — i — i — i — 5 — i — i — i — i — | — i — i — i — i — j — i — i — i — r 




I I I i I I I I i I L I I i I I I I i I I I I I 

50 100 150 200 850 

L (0-') 

Figure 2.1 Evolution of the kinetic, gravitational, and thermal energy per unit area, normalized to 
G 2 T,l/n 2 , in the standard run, which has L = 320GS o /fi 2 , resolution 1024 2 , and r co = 9.0x H^fT 1 . 

4 The natural unit that can be formed from G, £ and f2. 



32 



The turbulent state transports angular momentum outward via hydro dynamic and gravitational 
shear stresses. The dimensionless gravitational shear stress is 



1 



9rav ~ <f£ C 2> ./_T~4ttG 



dz 



9x9y 



(2.21) 



where g is the gravitational acceleration, and the dimensionless hydrodynamic shear stress is 



Qhyd 



(2.22) 



where () denote a spatial average. Figure 12.21 shows the evolution of (a grav ) and (cthyd) in the 
standard run. Averaged over the last 2301T 1 of the run, {{o-hyd)) = 0.0079, ((a grav )) = 0.017, and 
so the total dimensionless shear stress is {{a)) = 0.025, where (()) denote a space and time average. 




1 1 1 1 1 1 1 — I I i — I 1 — i i 

- gravitational 
hydrodynamic 



_i i i i i i_ 



... 







50 100 150 



200 



250 



Figure 2.2 Evolution of the gravitational and hydrodynamic pieces of (a) in the standard run. 
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The mean stability parameter (Q) = (c s )0/7rG ! (E) averages 1.86 over the last 230S7 1 of the 
run. Because the temperature and surface density vary strongly, other methods of averaging Q will 
give different results. 

Figure |231 sh ows a snap shot of the surface density at t = 50O 1 . The structure is similar to that 



observed in 



Gam mie (2001), with trailing density structures. The density structures are stretched 



into a trailing configuration by the prevailing shear flow. Their sc ale is determin ed by the disk 
temperature and surface density rather than the size of the box (see 



Gammic 



2001) 




Figure 2.3 Map of surface density at t = 50Q 1 in the standard run. Dark shades (blue in color 
version) indicate low density (0.2E o ) and light shades (yellow in color version) indicate high density 
(3S„). 



2.5.2 Varying E and ft 

We now turn to exploring the two-dimensional parameter space of models. First consider a series of 
models with the same initial central temperature, but with varying r co . As r co is lowered the time- 
averaged gravitational potential energy per unit area ((E g )) increases monotonically in magnitude. 
The gravito-turbulent state becomes more extreme, with larger {{a}}, larger perturbed velocities, 
and larger density contrasts. Eventually a threshold is crossed and the disk fragments. 
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Fragmentation is illustrated in Figure 12.41 which shows a snapshot from a run with Xl — 
6.6 x 10 3 gcm~ 2 , n = 5.4 x lCT^sec -1 . This corresponds to T Q = 1200, r co = 0.025ft -1 . The run 
has numerical resolution 256 2 and L = 80GT, o /Q 2 . The largest bound object in the center of the 
figure was formed from the collision and coalescence of several smaller bound objects. A snapshot 
of the optical depth at the same point in the simulation is given in Figure l2"31 For each snapshot, 
red indicates high values of the mapped variable and blue indicates low values. Much of the disk is 
optically thick, but most of the low density regions are optically thin in the Rosseland mean sense. 




Figure 2.4 Map of surface density in a run 
with T co = 0.025ft" 1 . Dark shades indicate 
both low density (10 _2 Xo, black in color ver- 
sion) and high density (10 2 £o, near the cen- 
ters of bound objects, red in color version). 




Figure 2.5 Map of optical depth r in a run 
with t co = 0.025ft -1 . Dark shades indicate 
both low t (10 -2 , black in color version) and 
high r (10 4 , near the centers of bound ob- 
jects, red in color version). 



Lowering r co sufficiently always leads to fragmentation. We have surveyed the parameter space 
of ft and S Q to determine where the disk begins to fragment. Each model was run to lOOft -1 . 5 
Figures 12.61 and 12.71 summarize the results. Two heavy solid lines are shown on each diagram. 
The upper line shows the most rapidly cooling simulations that show no signs of gravitational 
fragmentation (nonfragmentation point). Quantitatively, we define this as the point at which the 
time-averaged gravitational potential energy per unit area is equal to -3G 2 £3/ft 2 . 6 The lower line 
shows the most slowly cooling simulations to show definite fragmentation (fragmentation point). 
Quantitatively, we define this as the point at which the gravitational potential energy per unit area 



5 In four cases we had to run the simulation longer to get converged results. 

6 — 3 is the potential energy per unit area of a wave at the critical wavelength in a Q = 1 disk with 5E/E = ■s/S/it. 



No bound objects are observed throughout the duration of these runs. 
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Figure 2.6 Location of the critical curves as a function of initial volume density and temperature (in 
cgs units). Each contour line is an order of magnitude change in r co , solid/dotted lines indicating 
positive/negative integer values of log(r co ). 




Figure 2.7 Location of the critical curves as a function of initial surface density and rotation 
frequency (in cgs units). Each contour line is an order of magnitude change in r co , solid/dotted 
lines indicating positive/negative integer values of log(r co ). The gap in the center of the plot is due 
to the discontinuous jump in the value of fi. 
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is equal to — 300G 2 £g/O 2 at some point during the run. 7 Figure 121)1 shows the data in the p ,T 
plane, while Figure l2~Tl shows the results in the X ,f2 plane. Light contours are lines of constant 

Tco- 

The transition from persistent, gravito-turbulent outcomes to fragmentation is gradual and 
statistical in nature. Figure 12.81 shows the gravitational potential energy per unit area in the 
transition region for a series of runs with T a = 1200 K. The abscissa is labeled with the initial 
cooling time t co Q. There is a gradual, approximately logarithmic increase in the magnitude of 
{(Eg)) as t co decreases. Runs in this region exhibit the transient formation of small bound objects 
which might collapse if additional physics (e.g. the effects of MHD turbulence) were included in 
the model. Eventually —((E g )) begins to increase dramatically, and we define the transition point 
as the beginning of this steep increase in gravitational binding energy. 




Figure 2.8 Mean gravitational potential energy as a function of initial cooling time for a series of 
models with varying initial cooling time and T Q = 1200. 

7 These runs exhibit bound objects that persist for the duration of the run. 
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Figure 12.91 shows the run of t co Q for the fragmentation point, transition point, and nonfrag- 
mentation point as a function of T Q . It is surprising that a disk can begin to exhibit signs of 
gravitational collapse for r co O as large as 10 6 , a nd evade collapse for r co fi as small as 0.02. A 



naive application of the results of 



Gammiei ([200 ll ) would suggest that fragmentation should occur 



for t co Q < 3. Evidently this estimate can be off by orders of magnitude, with the largest error for 
T 10 3 K, just below the opacity gap. 



10" F 



~l 1 — I I I M 



~1 1 1 I I I I I 



1 1 1 — I I I I I I 



~\ 1 — i — r~n 




j i i 



Figure 2.9 Initial cooling times at the points of non-fragmentation, fragmentation and transition. 



he physical argument for fragmentation at short cooling times is as follows (e.g. IShlosman et al 
1990). Thermal energy is supplied to the disk via shocks. Strong shocks occur when dense clumps 
collide with one another; this occurs on a dynamical timescale ~ Q^ 1 . If the disk cools itself 
more rapidly then shock heating cannot match cooling and fragmentation results. This argument 



38 



is apparently contradicted by Figure l2~§l The resolution lies in finding an appropriate definition of 
cooling time. The disk loses thermal energy on the effective cooling timescale 



((Tc))- 1 = 



((A)) 
«E0>' 



(2.23) 



Figure . 1UI shows the run of ((r c )) at the fragmentation, transition, and non-fragmentation points. 
Evidently ((t c )) at transition lies between O -1 and WQ^ 1 . Figure 12.111 shows the run of r co and 
((t c )) on the transition line. Just below the opacity gap they differ by as much as four orders of 
magnitude. 





Figure 2.10 Effective cooling times at the Figure 2.11 Initial and effective cooling times 
points of non-fragmentation, fragmentation at the transition between non-fragmentation 
and transition. and fragmentation. 

Why do r co and ((r c )) differ by such a large factor? The answer is related to the existence of 
sharp variations in opacity with temperature. Consider a disk near the lower edge of the opacity 
gap. Once gravitational instability sets in, fluctuations in temperature move parts of the disk into 
the opacity gap. There, the opacity is reduced by orders of magnitude. Since the cooling rate for 
an optically thick disk is proportional to the cooling time drops by a similar factor. Relatively 
small var iations in temp erature can thus produce large variations in cooling rate. 



As in 



Gammie (2001 



vation implies that 



), the result {{t c ))£1 > 1 also implies a constraint on {(a)). Energy conser- 
3. 



n((w xy )) = ((A)), 



(2.24) 
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where W xy is the total shear stress (hydro dynamic plus gravitational). Equivalently, stress by 
rate-of-strain is equal to the dissipation rate. Using the definition of ((r c )), this implies 

««»= (7(7 -l)|n«r e ») • (2.25) 

Hence ((r c ))0 > 1 implies ((a)) < 1. Figure 12". 121 shows ((a)) vs ((r c )) for a large number of runs 
plotted against equation ()2.25j) . For small values of ((r c )) the numerical values lie below the line. 
These models are not in equilibrium (i.e., not in a statistically-steady gravito-turbulent state), so 
the time average used in equation ()2.24j) is not well defined. For larger values of ((t c )) numerical 
results typically (there is noise in the measurement of both ((a)) and ((r c )) because the time average 
is taken over a finite time interval) lie slightly above the analytic result. 




i 10 100 

«t c » (Q~ l ) 

Figure 2.12 Time-averaged shear stress vs. effective cooling time for a series of runs. The solid line 
shows the analytic result, based on energy conservation, from equation (|2.25[) . 
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The bias toward points lying slightly above the line reflects the fact that ((a)) measures the 
rate of energy extraction from the shear while ((r c )) measures the rate at which that energy is 
transformed into thermal energy. If energy is lost, perhaps to numerical averaging at the grid 
scale, then more energy must be extracted from the shear flow to make up the difference. Overall, 
however, the agreement with the analytic result is good and demonstrates good energy conservation 
in the code. 

The relationship between ((r c )) and ((a)) is interesting but not particularly useful because 
((r c )) is no more readily calculated than ((a)); it depends on a complicated moment of the surface 
density and temperature. Only for constant cooling time have we been able to evaluate this moment 
analytically. 

2.5.3 Isothermal Disks 

We have assumed that external illumination of the disk is negligible. This approximation is valid 
when the effective temperature Tj rr of the external irradiation is small compared to the central 
temperature of the disk. In the opposite limit, illumination controls the energetics of the disk and 
it is isothermal (if it is illuminate d directly so that shadowing effects, such as those considered by 



Jang-Condell and Sasselovl (|2003T ) are negligible). 

It is therefore worth studying the outcome of gravitational instability in an isothermal disk. 
The isothermal disk model has a single parameter: the initial value of Q. We ran models with 
varying values of Q and with (5v 2 )/c1 = 0.1. We find that models with Q < 1.4 fragment. 

It is likely that the mass of the fragments, etc., depends on how an isothermal disk becomes 
unstable. Rapid fluctuation of the external radiation field is likely to produce a different outcome 
than dimming on a timescale long compared to the dynamical time. 



2.6 Discussion 



Using numerical experiments, we have identified those disks that are likely to fragment absent 
external heating. Disks with effective cooling times ((t c )) < ar e susceptible to f ragme ntation. 



Shlosman et al 



1 199ft : if the 



This is what one might expect based on the simple argument of 
disk cools more quickly than the self-gravitating condensations can collide with one another, then 
those collisions (which occur on a timescale ~ cannot reheat the disk and fragmentation is 

inevitable. But our results are at the same time surprising. 
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The effective cooling time depends on the nonlinear outcome of gravitational instability. It 
depends on the cooling function, which in turn depends sensitively on S and U. Since £ and U 
vary strongly over the disk once gravitational instability has set in, it is difficult to estimate ((t c )) 
directly. One might be tempted to estimate ((r c ))(S, fl) ~ r C0 (S , Q, Q = 1), but our experiments 
show that this estimate can be off by as much as four orders of magnitude. The effect is particularly 
pronounced near sharp features in the opacity. For example, consider a model initially located just 
below the opacity gap with t co SI 3> 1. Gravitational instability creates dense regions with higher 
temperatures, where dust is destroyed. The result is rather like having to shed one's blanket on 
a cold winter morning: the disk loses its thermal energy suddenly. Pressure support is lost and 
gravitational collapse ensues. 

The difference between ((r c )) and r co (Q = 1) implies that a much larger region of the disk is 
susceptible to fragmentation than naive estimates based on the approximation ((r c )) ~ r co might 
suggest. For example, consider an equilibrium disk model with Q 3> 1 at small r. As r increases, Q 
declines. Eventually Q ~ 1 and gravitational instability sets in. There is then a range of radii where 
Q ~ 1, ((r c ))fi > 1 and recurrent gravitational instability can transport angular momentum and 
prevent collapse. Generally speaking, however, the cooling time decreases with increasing radius. 
Eventually ((r c ))f2 ~ 1 and fragmentation cannot be avoided. By lowering our estimate of ((t c )), 
we narrow the range of radii over which recurrent gravitational instability can occur. 

The general sense of our result is that it is extremely difficult to prevent a marginally-stable, 
Q ~ 1, optically-thick disk from fragmenting and forming planets (in circumstellar disks) or stars 
(in circumstellar and circumnuclear disks). This is particularly true for disks with T ~ 10 3 K, 
whose opacity is dominated by dust grains, i.e. disks whose temperature lies within a factor of 
several of the opacity gap. 

Our numerical model uses a number of approximations. First, our treatment is razor-thin, i.e. 
all the matter is in a thin slice at z = 0. The e f fect o f finite thickness on linear stability has 



been understood since 



Goldreich and Lvnden-Belll (1965): it is stabilizing because gravitational 



attraction of neighboring columns of disk is diluted by finite thickness. The size of the effect 
may be judged by the fact that Q = 0.676 is required for marginal stability of a finite-thickness, 
isothermal disk. 

The behavior of a finite-thickness disk in the nonlinear regime is more difficult to predict. Shocks 
will evidently deposit some of their energy away from the midplane, where it can be radiated away 
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more quickly (because the energy is deposited at smaller optical depth - see lPickett et alJ l2000) . 
Radiative diffusion parallel to the disk plane (not included here) may enhance cooling of dense, hot 
regions. Both these effects are destabilizing. Ultimately, however, a numerical study is required. 
This is numerically expensive: one must resolve the disk vertically, on the scale height H, and 
horizontally, at the critical wavelength 2nQH. 

Second, we have ignored magnetic fields. While there may be astrophysical situations where 
cool disks have such low ionization that they are unmagnetized, most disks are likely to contain 



dyna: 
(e-g- 



mically important magnetic fields that give rise to a dimensionless shear stress {{a)) > 0.01 



Hawlev et al 



1995). These fields are likely to remove spin angular momentum from partially 
collapsed objects, destabilizing them. Numerical experiments including both gravitational fields 



Balbu s and Hawlev 



and m agnetohydrodynamics are necessarily three dimensional (the instability of 
(1991) requires d z ^ 0), and are thus numerically expensive. 

Third, we have fixed 7 and \i for the duration of each simulation. This eliminates the soft spots 
in the equation of state associated with ionization of atomic hydrogen and dissociation of molecular 
hydrogen. In these locations the three dimensional 7 dips below 4/3, which is destabilizing. 

fourth, we h a ve tre ated the physics of grain destruction and formation very simply. In using 



the 



Bell and Lid (199J) opacities we implicitly assume that grains reform in cooling gas on much 



less t han a dynamical time. It is likely that grain re-formation will take some time (e.g. 



Hessman 



1991) and this will further reduce the disk opacity and enhance fragmentation. 



Fifth, we have neglected the effects of illumination. In the limit of strong external illumination, 
i.e. when the effective temperature of the irradiation Tj rr is large compared to the disk central 
temperature T c , the disk is isothermal (here T c is the temperature of a dense condensation). We 
have carried out isothermal experiments and shown that, for initial velocity perturbations with 
(5v 2 )/c 2 = 0.1, disks with Q < 1.4 fragment. Weaker illumination produces a more complicated 
situation that we have not explored here. Illumination-dominated disks that become unstable 
presumably do so because the external illumination declines, and the rate at which the external 
illumination changes may govern the nonlinear outcome. 

We conclude that disks with ((t c ))Q < 1 do not exist. Cooling in this case is so effective that 
fragmentation into condensed objects- stars, planets, or smaller accretion disks- is inevitable. 
As an example application of this result, consider the model for the nucleus of NGC 1068 

( 2003|) . Their model is an extended marginally-stable self- 



recently proposed by 



Lodato and Bertin 
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gravit ating disk of the type i nvestigated here a nd originally proposed by 



Goldreich and Lvnden-Belll 



( 19(55 ) for galactic disks and Pacz vnskl 1)1978 ) for accretion disks, although their disk is sufficiently 
massive that it modifies the rotation curve as well. Based on their Figure 3, at a typical radius 
of 0.5 pc, T, Q ~ 10 4 and £2 ~ 10~ 9 . According to our Figure 7 this disk is about 2 orders of 
magnitude too dense to avoid fragmentation. While it may be possible to avoid this conclusion 



by inv oking strong external h eating, the energy r e quirem ents are severe, as outlined in 



1 2003 ). The disk proposed by 
short timescale. 



Lodato and Bertin 



Goodman 



(2003) would therefore fragment into stars on a 
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Linear Theory of Thin, 
Radially-Stratified Disks 



3.1 Chapter Overview 

We consider the nonaxisymmetric linear theory of radially-stratified disks. We work in a shearing- 
sheet-like approximation, where the vertical structure of the disk is neglected, and develop equa- 
tions for the evolution of a plane-wave perturbation comoving with the shear flow (a shearing 
wave, or "shwave"). We calculate a complete solution set for compressive and incompressive short- 
wavelength perturbations in both the stratified and unstratified shearing-sheet models. We develop 
expressions for the late-time asymptotic evolution of an individual shwave as well as for the expec- 
tation value of the energy for an ensemble of shwaves that are initially distributed isotropically in 
fc-space. We find that: (i) incompressive, short-wavelength perturbations in the unstratified shear- 
ing sheet exhibit transient growth and asymptotic decay, but the energy of an ensemble of such 
shwaves is constant with time; (ii) short-wavelength compressive shwaves grow asymptotically in 
the unstratified shearing sheet, as does the energy of an ensemble of such shwaves; (iii) incompres- 
sive shwaves in the stratified shearing sheet have density and azimuthal velocity perturbations ST,, 
Sv y ~ t~ Rl (for |Ri| <C 1), where Ri = N^/{qVL) 2 is the Richardson number, N% is the square of the 
radial Brunt- Vaisala frequency and qfl is the effective shear rate; (iv) the energy of an ensemble 
of incompressive shwaves in the stratified shearing sheet behaves asymptotically as Rit 1-4Rl for 
|Ri| <C 1. For Keplerian disks with modest radial gradients, |Ri| is expected to be <C 1, and there 
will therefore be weak growth in a single shwave for Ri < and near-linear growth in the energy 

of an ensemble of shwaves, independent of the sign of Ri. 1 

lr To be published in ApJ Volume 626, Issue 2. Reproduction for this dissertation is authorized by the copyright 
holder. 
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3.2 Introduction 



Angular momentum transport is central to the evolution of astrophysical disks. In many disks 



angular momentum is likely redistributed internally b y magnetohydrodynamic 



driven by the magnetorotational instability (MRI; see 



Balbus and Hawlev 



MHD) turbulence 



1998). But in portions 



of disks around young, 



transients in quiescence ([Stone et al 



ow-mass stars, in cataclysmic- variable disks in quiescence , and in X-ray 



2000; 



Gammie and Me nou 



1998; 



Menou . 



2000), disks may be 



composed of gas that is so neutral that the MRI fails. It is therefore of interest to understand if there 
are purely hydrodynamic mechanisms for driving turbulence and angular momentum transport in 
disks. 

The case for hydrodynamic angular momentum transport is not promising. Numerical experi- 
ments carried out under conditions similar to those under which the MRI produces ample angular 
momentum fluxes- local shearing-bo x models- show small or n egative angular momentum fluxes 



when the magnetic field is turned off (H awlev et al 



1995 



1996). Unstratified shearing-sheet mod- 



els show decaying angular momentum flux and kinetic energy when nonlinearly perturbed, yet 
recover the well known, high Reynolds number nonlinear inst ability of plane Couette flow when the 



pa rameters of the mode 



by 



Umurhan and Regev 



are s et appropriately (B albus et al 



1996; see, however, the recent results 



2004). Local models with unstable vertical stratification show overturning 



and the developm ent of convective turbulen ce, but the mean angular momentum flux is small and 



of the wrong sign ( St one and B albu 



1996). 



Linear theory of global disk models has long indicat ed the presence of instabilitie s associated 
with reflecting b oundaries or features in the flow (see e.g. 



Gold reich et al 



1986; 



Goodman et al 



1987; 



Papaloizou 



Naravan et al.1 



1987 



and Pringle 



1984 



Lovelace et al 



1985, 



1999 



1987 



Li et al 



2000). Numerical simulations of the nonlinear outcome of th ese instabilit i es suggest tha t they 



1987 



Hawlev . 



199J). One 



saturate at low levels and are turned off by modest accretion (jBlaesJ, 
might guess that in the nonlinear outcome these instabilities will attempt to smooth out the fea- 
tures that give rise to them, much as convection tends to erase its parent entropy gradient. There 
are some suggestions, however, that such instabilities saturate into long-lived vort ices, which ma y 



2001). 



serve as obstructions in the flow that give rise to angular momentum transport (|Li et al 
We will consider this possibility in a later publication. 

Linear theory has yet to uncover a local instability of hydrodynamic disks that produces 
astrophysically-relevant angular momentum fluxes. Because of the absence of a complete set of 
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modes in the shearing-sheet model, however, local linear stability is difficult to prove. Local non- 
linear stability may be impossible to prove. Comparison with laboratory Couette flow experiments 
is complicated by several factors, not least of which is the inevitable presence of solid radial bound- 
aries in the laborator y that have no analogue in astro physical disks. 

Recently, however. iKlahr and Bodenheimerl (|2003) (hereafter KB03) have claimed to find a local 
hydrodynamic instability in global numerical simulations: the "Global Baroclinic Instability." The 
instability arises in a model with scale-free initial conditions (an equilibrium entropy profil e that 



varies as a power-law in radius) and thus does not depend on sharp features in the flow. 



Klahr 



(2004) has performed a local linear stability analysis of a radially-stratified accretion disk in an 
effort to explain the numerical results obtained by KB03. The instability mechanism invoked is the 
phenomenon of transient amplification as a shearing wave goes from leading to trailing. This is the 



mechanism that operates for nonaxisymme tric shearing waves in a dis 



c that i s nearly unstable to 



the axisymmetric gravitational 


1966 




Goldreich and Tre 


maine, 


linear analysis of 


KlahT 


(2004) 



1965; 



Julian and Toomre , 



1978). It is the purpose of this work to clarify and extend the 



low-ionization disks. 

To isolate the cause for instabilities originally observed in global 3D simulations, KB03 perform 
both local and global 2D calculations in the (R, </>)-plane. The local simulations use a new set 
of boundary conditions termed the shearing-disk boundary conditions. The model is designed to 
simulate a local portion of the disk without neglecting global effects such as curvature and horizontal 
flow gradients. The boundary conditions, which are described in more detail in KB03, require the 
assumption of a power-law scaling for the mean values of each of the variables, as well as the 
assumption that the fluctuations in each variable are proportional to their mean values. The radial 
velocity component in the inner and outer four grid cells is damped by 5% each time step in order 
to remove artificial radial oscillations produced by the model. 2 

The equilibrium profile for KB03's 2D runs was a constant surface density £ with either a 
constant temperature T or a temperature profile T oc R . The constant-T runs showed no 
instability while those with varying T (and thus varying entropy) sustained turbulence and positive 
Reynolds stresses. 3 The fiducial local simulations were run at a resolution of 64 2 , with a spatial 



2 It is not surprising that shearing disk boundary conditions as implemented in KB03 produce features on the 
radial boundary, because the Coriolis parameter is discontinuous across the radial boundary. 

3 Notice that with a constant S, the constant-T runs have no variation in any of the equilibrium variables, so it is 
not clear that the effects being observed in the 2D calculations are due to the presence of an entropy gradient rather 
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domain of R = 4 to 6 AU and A<fi = 30°. The unstable run was repeated at a resolution of 128 2 , 
along with a run at twice the physical size of the fiducial runs. One global model (with nonreflecting 
outflow boundary conditions) was run at a resolution of 128 2 with a spatial domain of R = 1 to 
10 AU and A<p = 360°. All the runs yielded similar results, with the larger simulations producing 
vortices and power on large scales. 

KB03 have chosen the term "baroclinic instability" by way of analogy with the baroclinic insta- 
bility that give s rise to weather patterns in the atmosphere of the Earth and other planets (see e.g. 



Pedloskv 



1979). 4 The analogy is somewhat misleading, however, since the baroclinic instability 



that arises in planetary contexts is due to a baroclinic equilibrium. In a planetary atmosphere, a 
baroclinically-unstable situation requires stratification in both the vertical and latitudinal direc- 
tions. 5 The stratification in KB03 is only in the radial direction, and as a result the equilibrium 
is barotropic. It is the perturbations that are baroclinic; i.e., the disk is only baroclinic at linear 



order in t 



Cabot 



re amp litude of a disturbance. 
and 



(198 



Knobloch and SpruitJ l[1986r ) have analyzed a thin disk with a baroclinic 



equilibrium state (with both vertical and radial gradients). The latter find that due to the dominant 
effect of the Keplerian shear, the instability only occurs if the radial scale height is comparable to 
the vertical scale height, a condition which is unlikely to be astrophysically relevant. As pointed out 
in KB03, the salient feature that is common to their analysis and the classical baroclinic instability 
is an equilibrium entropy gradient in the horizontal direction. As we show in §2, however, an 
entropy gradient is not required in order for two-dimensional perturbations to be baroclinic; any 
horizontal stratification will do. 

The "Global Baroclinic Instability" claimed by KB03 is thus analogous to the classical baroclinic 
instability in the sense that both have the potential to give rise to convection. 6 When neglecting 
vertical structure, however, the situation in an accretion disk is more closely analogous to a sh earing , 



stratified atmosp 



1961 



Chimonas 



iere, the stability of which is governed by the classical Richardson criterion (Miles, 



19701 ) . The only additional physics in a disk is the presence of the Coriolis force. 



Most analyses of a shearing, stratified atmosphere, however, only consider stratification profiles that 



than due simply to the presence of a pressure gradient. 

4 A baroclinic flow is one in which surfaces of constant density are inclined with respect to surfaces of constant 
pressure. If these surfaces c oincide , the fl ow is termed barotro pic. 

5 Contrary to the claim in lKlahil l)2004t) . the two-layer model l)PedloskvMl979f) does not ignore the vertical structure; 
it simply considers the lowest-order vertical mode. 

The classical baroclinic instability gives rise to a form of "sloping convection" llHoughtorl 12002ft since the latitu- 
dinal entropy gradient is inclined with respect to the vertical buoyancy force. 
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are stable to convection. The primary question that iKlahrl (J2004]) and this work are addressing, 
then, is whether or not the presence of shear stabilizes a stratified equilibrium that would be 
unstable in its absence. 

We begin in i!3.3l bv outlining the basic equations for a local model of a thin disk. § N3.4l and l3.5l 
describe the local linear theory for nonaxisymmetric sinusoidal perturbations in unstratified and 
radially-stratified disks, respectively. We summarize and discuss the implications of our findings in 

EH 

3.3 Basic Equations 

The effect of radial gradients on the local stability of a thin disk can be analyzed most simply 
in the two-dimensional shearing-sheet approximation 7 . This is obtained by a rigorous expansion 
of the equations of motion in the ratio of the vertical scale height H to the local radius R, fol- 
lowed by a vertical integration of the fluid equations. The basic equations that one obtains (e.g., 



Goldrei ch and Trem ainc 



IflZfiP are 

f + EV-.-O, (3.1) 



dv VP _ ^ 9 . . 

— + -— + 20 x v - 2qQ 2 xx = 0, (3.2) 

(XL z-j 



d InS 

— = 0, (3.3) 

where E and P are the two-dimensional density and pressure, S = -PE~ 7 is monotonically related 
to the fluid entropy, 8 v is the fluid velocity and d/dt is the Lagrangian derivative. The third 
and fourth terms in equation Q3.2j) represent the Coriolis and centrifugal forces in the local model 
expansion, where f2 is the local rotation frequency, x is the radial Cartesian coordinate and q is 
the shear parameter (equal to 1.5 for a disk with a Keplerian rotation profile). The gravitational 
potential of the central object is included as part of the centrifugal force term in the local-model 
expansion, and we ignore the self-gravity of the disk. 



7 See Rvu and Goo dman! Jl992f) fo r a dis cussion of why this approximation is appropriate for an analysis of local 
stability. See als"o lMaTcu^and^ressl il9T7f) . who use a similar approach to demonstrate the stability of unbounded 
viscous plane Couette flow. 

8 With the assumptions of vertical hydrostatic equilibrium and n egligible self-gravity, th e effective two-dimensional 
adiabatic index can be shown to be 7 = (37313 — 1)/(73D + 1) (e.g. iGoldreich et al.lll986t) . 
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It is worth emphasizing at this point that we have integrated out the vertical degrees of freedom 
in the model. We will later focus on perturbations with planar wavelengths that are small compared 
to a scale height, and these perturbations will be strongly influenced by the vertical structure of 
the disk. 

Equations ()3.1j) through (|3.3|) can be combined into a single equation governing the evolution 
of the potential vorticity: 

d (V x v + 2Sl\ _ d£ VExVP 



dt \ S J dt S 3 

In two dimensions, £ has only one nonzero component and can therefore be regarded as a scalar. 
Equation (|3.4[) demonstrates that for P = -P(E) (as in the case of a strictly adiabatic evolution with 
isentropic initial conditions), the potential vorticity of fluid elements is conserved. For P ^ -P(S), 
however, the potential vorticity evolves with time. A barotropic equilibrium stratification can result 
in baroclinic perturbations that cause the potential vorticity to evolve at linear order. This can be 
seen by linearizing the scalar version of equation (|3.4j) : 

— + v -V6{ + 6v- V£ = ^3 , (3-5) 

where we have dropped the term oc V£o X VPo- Notice that an entropy gradient is not required 
for the evolution of the perturbed potential vorticity. For So = -Po^o = constant, equation (J3,5|) 
reduces to 

-^ + v -V5t + 6vVZ = — V " '-. 3.6 

at 7 s ipo 

Potential vorticity is conserved only in the limit of zero stratification (Pq = constant) or adiabatic 
perturbations (5S = 0). 

3.4 Unstratified Shearing Sheet 

Our goal is to understand the effects of radial stratification, but we begin by developing the linear 
theory of the standard (unstratified) shearing sheet, in which the equilibrium density and pressure 
are assumed to be spatially constant. This will serve to establish notation and method of analysis 
and to highlight the changes i ntroduced by radial stratificatio n in the next section. 



Our analysis follows that of 



Goldreich and Tremainel (jl97Sl ) except for our neglect of self-gravity. 
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The equilibrium consists of a uniform sheet with S = So = constant, P = Pq = constant, and 
vq = —q£lxy. We consider nonaxisymmetric Eulerian perturbations about this equilibrium with 
space-time dependence 5(t)exp(ik x (t)x + ik y y), where 

k x (t) = k x0 + qilkyt (3.7) 

(with k x o and k y > constant) is required to allow for a spatial Fourier decomposition of the 
perturbation. We will refer to these perturbations as shearing waves, or with some trepidation, but 
more compactly, as "shwaves". 

3.4.1 Linearized Equations 

To linear order in the perturbation amplitudes, the dynamical equations reduce to 

— — h ik x 5v x + ikySvy = 0, (3-8) 
So 



5P 

5v x - 2fl5v y + ik x — = 0, (3.9) 
So 



5P 

5v y + (2 - q)Sl8v x + ik y — = 0, (3.10) 

So 



5P 

— + c 2 s (ik x 5v x + ikySvy) = 0, (3.11) 

where c 2 s = 7P0/S0 is the square of the equilibrium sound speed and an over-dot denotes a time 
derivative. 

The above system of equations admits four linearly-independent solutions. Two of these are 
the nonvortical shwaves (solutions for which the perturbed potential vorticity is zero), which in 
the absence of self-gravity can be solved for exactly. The remaining two solutions are the vortical 
shwaves. When k y — ► the latter reduce to the zero-frequency modes of the axisymmetric version 
of equations (|3.8|) through (|3.11|) . One of these (the entropy mode) remains unchanged in nonax- 
isymmetry (in a frame comoving with the shear). There is thus only one nontrivial vortical shwave 
in the unstratified shearing sheet. 

In the limit of tightly- wound shwaves {\k x \ 3> k y ), the nonvortical and vortical shwaves are 
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compressive and incompressive, respectively. In the short-wavelength limit {Hk y 3> 1, where 
H = c s /Vt is the vertical scale height), the compressive and incompressive solutions remain well 
separated at all times, but for Hk y < 0(1) there is mixing between them ne ar k x = as an 



1997, ; also Goodman 



incompressive shwave shears from leading to trailing. (Chagelish vili et al. 
2005, private communication) With the understanding that the distinction between compressive 
shwaves and incompressive shwaves as separate solutions is not valid for all time when Hk y < 
0(1), we generally choose to employ these terms over the more general but less intuitive terms 
"nonvortical" and "vortical." 

Based upon the above considerations, it is convenient to study the vortical shwave in the short- 
wavelength, low-frequency (dt <C c s k y ) limit. This is equivalent to working in the Boussinesq 
approximation, 9 which in the unstratified shearing sheet amounts to assuming incompressible flow. 
In this limit, equation 1)3.8)1 is replaced with 

k x 5v x + kySvy = 0. (3-12) 

This demonstrates the incompressive nature of the vortical shwave in the short- wavelength limit. 

3.4.2 Solutions 

In the unstratified shearing sheet, equation ()3.5)) for the perturbed potential vorticity can be inte- 
grated to give: 

ik x 5v y - ik y 5v x 5Y, 
o€u = ^ £ot^ = constant, (3.13) 

where £o = (2 — q)£l/T>o is the equilibrium potential vorticity and we have employed the subscript 
u to highlight the fact that the perturbed potential vorticity is only constant in the unstratified 
shearing sheet. To obtain the compressive-shwave solutions, we set the constant <5£ n to zero. 
Combining equations 1)3. lOj) and 1)3.13)1 with 5t; u = 0, one obtains an expression for 5v xc in terms 
of 5v yc and its derivative: 

* _ c 2 s k x k y 5v yc - CogoHc r „ 1A x 

0Vxc ~ t2y2 , „2J.2 ' 

where the subscript c indicates a compressive shwave. The associated density and pressure pertur- 
bations are 

,„ _ SP C _ Xo^okxSvyc + k y 5v yc 

6 ^ C -~^T- 1 p2y2 T C 2U2 \y- li3 ) 



9 We demonstrate this equivalence in the Appendix. 
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via equation (|3.13|) . Reinserting equation (|3.14|) into equation ()3.1(Jj) . taking one time derivative 
and replacing 5P via equation ()3.11() . we obtain the following remarkably simple equation: 



Sv yc + (c 2 s k 2 + k 2 ) 5v yc = 0, 



(3.16) 



where k 2 = k 2 + k 2 and k 2 = (2 — q)£l 2 is the epicyclic frequency. Changing to the dimensionless 
dependent variable 

h.r.h... f h.r.k... 

(3.17) 



and defining 



rp _ • j^Csky I k X Q 



2C^ A^-y 

«4/ 



(3.18) 



the equation governing 5v yc becomes 



d 5v yc ( 1 



dT 2 



+ [-T z -C) 5v yc = 0. 



(3.19) 



This is the parabolic cylinder equation (e.g. lAbramowitz and Ste gun 1972), the solutions of which 
are parabolic cylinder functions. One representation of the general solution is 



Sv 



]_rp2 

e 2 



ClM [--^C,-,-T 2 ) +C2 Tm(--^C,-, 1 -T 2 
l 4 2 ' 2' 2 / V 4 2' 2' 2 



(3.20) 



where c\ and C2 are constants of integration and M is a confluent hyper geometric function. This 
completely specifies the compressive solutions for the unstratified shearing sheet, for any value of 



k y . 



Equation (|3.19|) has been analyzed in detail by 



Naravan et al 



( 19871 ): their modal analysis yields 



the analogue of equation (|3.19|) in radial-position space rather than in the radial-wavenumber 
(k x = k y r) space that forms the natural basis for our shwave analysis. One way of seeing the 
correspondence between the modes and shwaves is to take the Fourier transform of the asymptotic 
form of the solution. Appropriate linear combinations of the solutions given in equation (|3.2()|) have 
the following asymptotic time dependence for r> 1: 



bv yc oc \\ — exp ( ±-T ) oc / _ exp ( ±i / c s k x dt 



1 

k x 



(3.21) 



The Fourier transform of the above expression, evaluated by the method of stationary phase for 
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Hky » 1, yields 



5v yc {X) oc 



(3.22) 



Naravan et al 



(1983), m 



which is equivalent to the expressions given for the modes analyzed by 
which the dimensionless spatial variable (with zero frequency, so that corotation is at x = 0) is 
defined as 

(3.23) 



X 



-x. 



To obtain the incompressive shwave, we use the condition of incompressibility (equation Q3.12[l 1 
to write 5v y in terms of 5v x , and then combine the dynamical equations l)3.9j) and (|3,1U|) to eliminate 
5P. The incompressive shwave is given by: 



k 2 

5v xi = 5v xi0 -j^, 



(3.24) 



Sv 



k :r 



Sv a 



(3.25) 



5Ei 
So 



6^ 
7P0 



ic s hy 



+ 2 (q-l)n — 



(3.26) 



where the subscript i indicates an incompressive shwave, k 2 



k 2 n + k 2 and 5v x io is the value of 



8v x i at t = 0. This solution is uniformly valid for all time to leading order in (Hk y ) <C 1. 



3.4.3 Energetics of the Incompressive Shwaves 

We define the kinetic energy in a single incompressive shwave as 



Eki = -Z (6v xi + 8v. 



V 2 



k 4 

V X 2 "'O 



(3.27) 



which peaks at k x = 0. This is not the o nly possible definition for the energy associated with a 
shear-flow disturbance; see Appendix A of 



Naravan et al 



(1987) for a discussion of the subtleties 



involved in defining a perturbation energy in a differentially-rotating system. The energy defined 
above can simply be regarded as a convenient scalar measure of the shwave amplitude. 



1( lChagelishvili et alJ |2003) obtained this solution by starting with the assumption of incompressibility. In the 
incompressible limit, it is an exact nonlinear solution to the fluid equations. 
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One can also define an amplification factor for an individual shwave, 



A 



E kl {k x = 0) 
Ekiit = 0) 



1 + 



k 2 . 



xO 

Ky 



(3.28) 



which indicates that an arbitrary amount of transient amplification in kinetic energy can be obtained 
as one increases the a mount of swing for a lead i ng shwave (k x n <C — hX T his is essentially the 



mecha nism invoked by 



C hagelishvili et al 



Umurhan and Regevl (j20oJ) and 



Afshordi et al 



(2003J), 

(200^) to argue for the onset of turbulence in unmagnetized Keplerian disks. 

Because only a small subset of all Fourier components achieve large amplification (those with 
initial wavevector very nearly aligned with the radius vector), one must ask what amplification 
is achieved for an astrophysically relevant set of initial conditions containing a superposition of 
Fourier components. It is natural to draw such a set of Fourier components from a distribution 
that is isotropic, or nearly so, when ko is large. 

Consider, then, perturbing a disk with a random set of incompressive perturbations (initial 
velocities perpendicular to fco) drawn from an isotropic, Gaussian random field and asking how the 
expectation value for the kinetic energy associated with the perturbations evolves with time. The 



evolution of the expected energy density is given by the following integral: 



(Ei) = L 2 j d 2 k (E ki ) = L 2 j d 2 k ^ (5v 2 xl0 ) 



k 1 
k 2 y k 2 ' 



(3.29) 



where () indicates an average over an ensemble of initial conditions, the first equality follows from 
Parseval's theorem, the second equality follows from the incompressive shwave solution (|3,24j) - (|3.26|) 
and therefore applies only for k^H 3> 1, and L 2 is a normalizing factor with units of length squared. 

For initial conditions that are isotropic in fco (<5^zio = 5v±(ko, 9) sin 9, where (5f^(A;o)) is the 
expectation value for the initial incompressive perturbation as a function of ko and tan 9 = k y /k x o), 
the integral becomes 



(Ek) = \z L 2 J k dk (8vl(k )) J\e 



1 



sin 9 + (qflt sin 9 + cos I 



(3.30) 



Changing integration variables to r = q£lt + cot 9, the angular integral becomes 



oo 2 

dr K - 2z 

-oo 



1 + T 2 



(3.31) 



55 



which is independent of time; hence 



(Ei) = {Ei(t = 0)) (3.32) 
and we do not expect the total energy in incompressive shw aves to evo l ve. 11 This same calculation 



Shepherd 



(1985), who also points out 



has been performed in the context of plane Couette flow by 
that the amplification factor due to a distribution of wavevectors in an angular wedge AO has an 
upper bound of 2tt/(A9). This indicates that the amplification will be modest unless the initial 
disturbance is narrowly concentrated around a single wavevector. 

Although this result may appear to depend in detail on the assumption of isotropy, one can 
show that it really only depends on (E}-i(t = 0)) being smooth near sine? = 0, i.e. that there should 
not be a concentration of power in nearly radial wavevectors. This can be seen from the following 
argument. If we relax the assumption of isotropy, the angular integral becomes 



I 



*« , <M^» 5 . (3.33) 

sin 2 (9 + (^ sin (9 + cos 0) 2 



For q£lt 1 the above integrand is sharply peaked in the narrow regions around tan 9 = —l/(qttt) <C 
1 (i.e., s'mO ~ 0). One can perform a Taylor-series expansion of (5v^(ko,9)} in these regions, and 
as long as (5v±(ko,0)) itself is not sharply peaked it is well approximated as a constant. A modest 
relaxation of the assumption of isotropy, then, will result in an asymptotically constant value for 
the energy integral. 

Based upon this analysis, large amplification in an individual shwave does not in itself argue for 
a transition to turbulence due to transient growth. One must also demonstrate that a "natural" set 
of perturbations can extract energy from the background shear flow. In the case of the unstratified 
shearing sheet, the energy of a random se t of incompressive perturba tions remains constant with 



time. This is consistent with the results of 



Umurhan and Regevl 1)2004 ) . who see asymptotic decay 



in linear theory, because they work with a finite set of wavevectors, each of which must decay 
asymptotically. 



11 Notice that while the energy of each individual shwave decays asymptotically, the energy of an ensemble does 
not. This is due to the spread of amplification factors in the spectrum of shwaves; some are amplified by very large 
factors while others are amplified very little. 
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3.4.4 Energetics of the Compressive Shwaves 

Here we calculate the energy evolution of the compressive shwaves for comparison purposes. We will 
consider the evolution of short-wavelength compressive shwaves in which only the initial velocity 
is perturbed, both for simplicity and for consistency with our calculation of the short-wavelength 
incompressive shwaves. As before, we will assume that the initial kinetic energy is distributed 
isotropically. 

We use the WKB solutions to equation QM.lfiJ) with Hk y ^> 1. With the initial density 
perturbation set to zero (consistent with our assumption of only initial velocity perturbations), the 
uniformly- valid asymptotic solution to leading order in (Hky) -1 is given by 

Sv yc = 5v yc0 J ^ cos(^ - Wq), (3.34) 



k 

Sv xc = -j^Svyc, (3.35) 



1 

5T, C = -jrrSvyc, (3.36) 

Cg Ky 

where the WKB eikonal is given by 

J Cskdt = ^hL J ^/TT^dT = ^ ( r ^/rh72 + ln(r+ V / l+^)) , (3.37) 



W 



with Wo being the value of W at t = 0. 13 

Using equation ()3.35|) . the energy integral for the compressive shwaves in the short- wavelength 
limit is 

(E c ) = L 2 J d 2 k (E kc ) = L 2 J d 2 k ^o(Sv 2 yc )^. (3.38) 



With initial velocities now parallel to fco (and again isotropic), this becomes 

(E c ) = \^Q l2 J kodko(5v\(ko)) d9 ^/sin 2 6 + (qttt sin 6 + cos#) 2 cos 2 (W - W ). (3.39) 



12 These solutions are the short- wavelength, /wg/i-frequency (dt ~ 0(c s k y )) limit of the full set of linear equations 
in the shearing sheet; see the Appendix. 

13 This is not the same WKB solution that is calculated in the tight-winding approximation by 
iGoldreich and Tremainel J1978T) : in that case c B k y /n 1, the opposite limit to that which we are considering here. 
The two WKB solutions match for r ^> 1 in the absence of self-gravity. We have verified the accuracy of this solution 
by comparing it to the exact solution with acceptable results, and it is valid to leading order for all time. 
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For qQt 3> 1, the angular integral is approximated by 



f d6 \ sin 6\ (1 + cos(2W -2W )) ~2qnt + x ^^cos(c s k qnt 2 -ir/ 4), (3.40) 
Jo V c s ko 

where the second approximation comes from employing the method of stationary phase. 14 In the 
short-wavelength limit, then, 

(E c (qQt > 1)) = 2qQt (E c (t = 0)). (3.41) 

Thus the kinetic energy of an initially isotropic distribution of compressive shwaves grows, presum- 
ably at the expense of the background shear flow. 

The fate of a single compressive shwave is to steepen into a weak shock train and then decay. 
The fate of the field of weak shocks generated by an ensemble of compressive shwaves is less clear, 
but the mere presence of weak shocks does not indicate a transition to turbulence. 



3.5 Radially- St ratified Shearing Sheet 

We now generalize our analysis to include the possibility that the background density and pressure 
varies with x; this stratification is required for the manifestation of a convective instability. In order 
to use the shwave formalism we must assume that the background varies on a scale L ~ H <C R so 
that the local model expansion (e.g., the neglect of curvature terms in the equations of motion) is 
still valid. 

With this assumption the equilibrium condition becomes 



where a prime denotes an x-derivative. One can regard the background flow as providing an effective 
shear rate 

qfl = -v' (3.43) 

that varies with x, in which case vq = — J x q(s)ds £ly. 

Localized on this background flow we will consider a shearing wave with k y L 3> 1. That is, we 

14 The first approximation breaks down near sin# = 0, but the contribution of these regions to the integral is 
negligible for qQ.t 3> 1, in contrast to the situation for incompressive shwaves. 
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will consider nonaxisymmetric short- wavelength Eulerian perturbations with spacetime dependence 
6(t) exp(i J k x (t, s)ds + ik y y + ik z z), where k y and k z are constants and 



k x (t,x) = k x0 + q{x)Q,k y t. 



(3.44) 



It may not be immediately obviou s that this is a valid expansion since the shwaves sit on top 



of a radially-varying background (see 



Toomre 



1969 for a discussion of waves in a slowly-varying 



background). But this is an ordinary WKB expansion in disguise. To see this, one need only 
transform to "comoving" coordinates x' = x, y' = y + f x q(s)ds£lt, t' = t (this procedure may be 



more familiar in a cosmological context; as 



Balbua (|1988n has pointed out, this is possible for any 



flow in which the velocities depend linearly on the spatial coordinates). In this frame the time- 
dependent wavevector given above is transformed to a time-independent wavevector. The price paid 
for this is that d x — > d x > + qVltd y i, so new explicit time dependences appear on the right hand side 
of the perturbed equations of motion, and the perturbed variables no longer have time dependence 
exp(iu;i'). Instead, we must solve an ODE for 5(t'). The y' dependence can be decomposed as 
exp(ik y y'). The x' dependence can be treated via WKB, since the perturbation may be assumed 
to have the form W(ex', et') exp(ik' ■ x'). This "nearly diagonalizes" the operator d x i. Thus we are 
considering the evolution of a wavepacket in comoving coordinates — a "shwavepacket" . 

For this procedure to be valid two conditions must be met. First the usual WKB condition 
must apply, k y L S> 1. Second, the parameters of the flow that are "seen" by the shwavepacket must 
change little on the characteristic timescale for variation of 5(t), which is O -1 for the incompressive 
shwaves. For solid body rotation (q = 0) the group velocity (derivable from equation |3.59j . below) 
is \v g \ < N x /k (for positive squared Brunt- Vaisala frequency N x , defined below; for N x < 
the waves grow in place), so the timescale for change of wave packet parameters in this case is 
L/\v g \ > kL/N x J2 . It seems reasonable to anticipate similarly long timescales when shear 
is present. As a final che ck, we have verified directly, using a code based on the ZEUS code of 
Stone and Normanl ()1992l ). that a vortical shwavepacket in the stratified shearing sheet remains 
localized as it swings from leading to trailing. 
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3.5.1 Linearized Equations 

To linear order in the perturbation amplitudes, the dynamical equations reduce to 



|S ^ 6v^_ + ^g Vx + ij^fo _|_ ik z 5v z = 0, 



(3.45) 



5P (? <5£ 

5v x - 2fl5v y + ik x — -p-— = 0, 

2^0 Z^o 



(3.46) 



5P 

5v y + (2 - q)Q5v x + ik y — = 0, 



(3.47) 



SP 

5v z + ik z — = 0, 



(3.48) 



where 



5P 



,5E 



,Sv 7 



2,0 ^0 J^s 



1 
Lp~ 



P' 

7^0 



1 1 

is is 



y" <?' 
S 75 



(3.49) 



(3.50) 



define the equilibrium pressure, density and entropy length scales in the radial direction. We have 
included the vertical component of the velocity in order to make contact with an axisymmetric 
convective instability that is present in two dimensions, after which we will set k z to zero. 

We will be mainly interested in the incompressive shwaves because the short-wavelength com- 
pressive shwaves are unchanged at leading order by stratification. We will therefore work solely in 
the Boussinesq approximation. 15 In addition to the assumption of incompressibility, this approx- 
imation considers 5P to be negligible in the entropy equation; pressure changes are determined 
by whatever is required to maintain nearly incompressible flow. The original Boussinesq approx 



imation ap plies only to incomp r essib 



(1930) and 



Spiegel and Veroni; 



Jeffreys 



e fluids. It was extended to compressible fluids by 

(1960). We show in the Appendix that it is formally equivalent 

to taking the short- wavelength, low-frequency limit of the full set of linear equations. From this 

viewpoint, assuming that Hk y 8P / Pq is of the same order as the other terms in the dynamical 

equations implies that 5P/Pq ~ (Hk y )~ l 5Y l /Y<Q, thus justifying its neglect in the entropy equation. 
15 We also drop the subscripts distinguishing between the compressive and incompressive shwaves. 
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We therefore replace equations Q3.45|) and ()3.49j) with 



k x 5v x + k y 5v y + k z 5v z = (3.51) 

and 

— - -A = 0. (3.52) 

Using equations (j3.47j) and (|3.52j) and the time derivative of equation (|3.51|) , one can express 
5v y and 8P in terms of 5v x and 5v x : 

6P_ = Jk x 5v x + 2(q- l)Qk y 5v x ^ 
Ylfl ky -\- k 2 

(-qk 2 + (q - 2)k 2 )Q5v x - k x k y 5v x 
Sv v = " I2-T2 • ( 3 - 54 ) 

Eliminating 5P in equation (|3.46j) via equation ([3.53ft gives 

k 2 5v x + 2(q - l)nk x k y 5v x = (JfeJ + Jfe|)(2fitfu w + (c^/L F )5S/S ), (3.55) 

where A; 2 = A; 2 + A 2 + fc 2 . Taking the time derivative of this equation and eliminating 5T, and 8v y 
via equations (|3.52|) and (|3.54l) . we obtain the following differential equation for 8v x : 

k 2 5v x + 4qnk x ky5v x + [k 2 y (N 2 + 2q 2 n 2 ) + k 2 (N 2 + R 2 )] 5v x = 0, (3.56) 

where h 2 = 2(2 — q)£l 2 is the square of the effective epicyclic frequency and 



N 2 ee (3.57) 



2 



is the square of the Brunt- Vaisala frequency in the radial direction 



n> 



Notice that N£, q and k are all functions of x and vary on a scale L ~ H. 
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3.5.2 Comparison with Known Results 

Setting k y = in equation ()3.56|) yields the axisymmetric modes with the following dispersion 



relation (for S(t) oc e lujt ): 

.2 l A t2 , ~2 



;,^2(^ + - 2 )- (3-58) 

This is the origin of the H0iland stability criterion: the axisymmetric modes are stable for N 2 + k 2 > 
0. In the absence of rotation this reduces to the Schwarzschild stability criterion: N 2 > is the 
necessary condition for stability. The effect of rotation is strongly stabilizing: if N 2 < —k 2 , as 
required for instability, then LgLp ~ H 2 ; pressure and entropy must vary on radial scales of order 
the scale height for the disk to be H0iland unstable. 

Notice that effective epicyclic frequency k 2 only stabilizes modes with nonzero k z . The stability 
of nonaxisymmetric shwaves with k z = (as in the mid-plane of a thin disk) is the open question 
that this work is addressing. In this limit and in the absence of shear the Schwarzschild stability 
criterion is again recovered: with k z = and q = in equation (J3.56|) the dispersion relation 
becomes 

^ 2 = J^n-^l (3-59) 

K x0 K y 

If there is a region of the disk where the effective shear is zero, a WKB normal-mode analysis will 
yield the above dispersion relation and there will be convective instability for N 2 < 0. It appears 
from equation (|3.56|) that differential rotation provides a stabilizing influence for nonaxisymmetric 
shwaves just as rotation does for the axisymmetric modes. Things are not as simple in nonax- 
isymmetry, however. The time dependence is no longer exponential, nor is it the same for all the 
perturbation variables. There is no clear cutoff between exponential and oscillatory behavior, so 
the question of flow stability becomes more subtle. 

As discussed in the introduction, the Boussinesq system of equations in the shearing-sheet model 
of a radially-stratified disk bear a close resemblance to the system of equations employed in analyses 
of a shearing, stratified atmosphere. A sufficient condition for stability in the latter case is that 

N 2 1 

Ri = 7% > T (3-60) 

K) 2 4 

everywhere in the flow, where Ri is the Richardson number, a measure of the relativ e im portance 



of buo yancy and shear. This stability criterion was originally proved by 



Miles! (|l96l|) and 



Howar 



( 196lj ) for incompressible fluids, and its extension to compressible fluids was demonstrated by 

62 



Chimonaa (|197Q). The stability criterion is based on a normal-mode analysis with rigid boundary 
conditions. Other than differences in notation (e.g., our radial coordinate corresponds to the vertical 
coordinate in a stratified atmosphere), the key differences in our system are: (i) the equilibrium 
pressure gradient in a disk is balanced by centrifugal forces rather than by gravity; (ii) the disk 
equations contain Coriolis force terms; (iii) most atmospheric analyses only consider an equilibrium 
that is convectively stable, whereas we are interested in an unstable stratification; (iv) we do not 
employ boundary conditions in our analytic model since we are only interested in the possibility of 
a local instability. 

The lack of boundary conditions in our model makes the applicability of the standard Richardson 
stability criterion in determining local stability somewhat dubious, since the lack of boundary 
conditions precludes the decomposition of linear disturbances into normal modes. The natural 
procedure for performing a local linear analysis in disks is to decompose the perturbations into 
sh waves, as we ha v e don e . 



Eliassen et al 



(1953) consider both stable and unstable atmospheres and analyze an initial- 
value problem by decomposing the perturbations in time via Laplace transforms. For flow between 
two parallel walls, they find that an arbitrary initial disturbance behaves asymptotically as i( a_1 )/ 2 
for -3/4 < Ri < 1/4, where 

a = Vl -4Ri, (3.61) 

which grows algebraically for Ri < 0. The disturbance grows exponentially only for Ri < —3/4. 
For a semi-infinite flow, the power-law behavior in time holds for —2 < Ri < 1/4, with exponential 
growth for Ri < —2. These results illustrate the importance of boundary conditions in determining 
stability. 

In the k z = limit that we are concerned with here, the correspondence between the disk and 
atmospheric models turns out to be exact in the shwave formalism. This is because the Coriolis force 
only appears in equation (|3.56|) via k 2 , which disappears when k z = 0. The equation describing the 
time evolution of shwaves in both a radially-stratified disk and a shearing, stratified atmosphere is 
thus 

k 2 5v x + 4qnk x k y 5v x + k 2 (N 2 + 2q 2 n 2 ) 5v x = 0. (3.62) 

We analyze the solutions to this equation in the following section. 17 

17 This equation is also obtained in a shwave analysis of interchange instability in a disk with a poloidal magnetic 
field f S pruit et, all llflQfJ l , with N% replaced by a magnetic buoyancy frequency. 
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3.5.3 Solutions 

Changing time variables in equation (|3.62|) to f = k x /k y , the differential equation governing 5v x 
becomes 

9x d 2 5v T d5v r .„ , 
(l + f 2 )-7^ + 4f— - + {m + 2)5v x = 0. (3.63) 

«T Z (IT 

The solutions to equation (|3.63|) are hyper geometric functions. With the change of variables z = 
— f 2 , equation (j3.63j) becomes 

dHv^ l-5z d5v x Ri + 2 
,(1 - z)^- + — - —Sv, = 0. (3.64) 



The hypergeometric equation ()Abramowitz and SteeunL 11972) 



. d 2 5v r , , . . d5v r , , 

z(l - z)—^- + [c-(a + b + l)z] —r^ - ab5v x = 3.65 
dz z dz 

has as its two linearly independent solutions F(a, b; c; z) and z 1 ~ c F(a — c + 1, b — c + 1; 2 — c; 2). 
Comparison of equations (|3.64|) and (|3.65|) shows that a = (3 — a)/4, 6 = (3 + a)/4 and c = 1/2, 
where a is defined in equation (|3.61|) . 

The general solution for Sv x is thus given by 

3 — a 3 + a 1 _ 2 \ n ~ & { ^ ~ a 5 + a 3 _ 2 



ft* = F ^— , — ; -; -f^j + C 2 f F , _; -; -f- j , (3.66) 

where C\ and C2 are constants of integration representing the two degrees of freedom in our 
reduced system. These two degrees of freedom can be represented physically by the initial velocity 
and displacement of a perturbed fluid particle in the radial direction. The radial Lagrangian 
displacement £ x is obtained from equation (|3.66|) by direct integration, 18 

C2 „ fl — a 1 + a 1 _ 2 \ Ci-^/3 — a 3 + a 3 _ 2 



£ x = Sv x dt = — — — F ,— — ;-;-f 2 )+^fF[ ,— — . (3.67) 

s J x qnm V 4 4 2 J V 4 4 2 / 

The solutions for the other perturbation variables can be obtained from equations (|3.51|) . (|3.52|) 
and (TT331 with k z = 0: 

Sv y = —fSv x , (3.68) 

18 In our notation, a subscript x or y on the symbol £ indicates a Lagrangian displacement, not a component of the 
potential vorticity, which is a scalar. 



64 



and 



5P 



Pq ic s ky 



So Lg 



d ( 5v x \ „,„ „ . foa- 

9T— — ] +2(g-l)— ^ 



(3.69) 



(3.70) 



(if \ c s 

It can be seen from the latter equation and the solution for Sv x that 5P/Pq remains small compared 
to 5v x /c s in the short- wavelength limit. This demonstrates the consistency of the Boussinesq 
approximation. 



The hypergeometric functions can be transformed to a form valid for large f (see 
I972I equations 15.3.7 and 15.1.1). An equivalent form of the solution for |f| ^ 1 is 



Abramowitz and Stegun 



Sv x = (C 1 V l + sgn(T)C 2 V 2 )\T\ 2 

2+3 _ / 3 + a 5 + a 



a 1 



+ 



(C1F3 + sgn(f )C 2 V A ) \f\~-F ( ! + -;-_), 



2 -■ 



where sgn(f) is the arithmetic sign of f and the constants Vi are given by 

r(i)r(f) 



rmr(-^) 



, v 2 



^3 = 



r(i)r 



r(|)r(f) 
r(§)r(-f) 



r(^)r(-±±^) ' '*-r(^a)r(^)' 

Expanding the above form of the solution for |f| 3> 1, we obtain 



(3.71) 



(3.72) 



Sv x = (dV! + sgn(f )C 2 V 2 ) |f|^ + (C1F3 + sgn(f)C 2 y 4 ) |^]~" + 0(f- 2 ). 



(3.73) 



An equivalent form of £ x for \f \ ~S> 1 is 



C 2 X\ C\X 2 
+ sgn(T) 



C2X3 



\f\~F 



3 — a 1 — a 



1 



a 1 
2'~?2 



+ 



. . CVXkN . , a+i _ / 3 + a 1 + a a 1 

^Ri +sgn(r) ^rJ |T| " 2 — ;1+ 2 ; -^ 



(3.74) 
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where the constants X, are given by 



v _ r(i)r(f) r(|)r(f) 



r( i±a )r( i±a) > < ~ r (aja) r (»ja) ' 

v - r (l)r(-f) Y _ r(|)r(-f) 
r (i^) r (i^) ' A4 " r (s^s) r (s^) ' 

Expanding ^ for |f| S> 1 yields 

^• = r^ffii +sgn(r) ^rj |r| + v ^mr +sgn(r) ^rj |r| +0(r } - (3 ' 76) 

The dominant contribution for each perturbation variable at late times is thus 

8P at 5v x ~ t 3 ^ , (3.77) 



a— 1 



5Soc^x~t" 3 ~ , (3.78) 
and 



0—1 



^ cx tfe-c ~ i 2 . (3.79) 

This leads to one of our main conclusions: the density and y-velocity perturbations will grow 
asymptotically for a > 1, i.e. Ri cx N x < 0. 19 For small Richardson number, however (as is 
expected for a Keplerian disk with modest radial gradients), a ~ 1 — 2Ri and the asymptotic 
growth is extremely slow: 

<5E ~ 5v y ~ t~ Ri . (3.80) 

In the stratified shearing sheet, the right-hand side of equation 1)3. 5 J) governing the evolution of 
the perturbed potential vorticity is no longer zero. The form of this equation for the incompressive 
shwaves is 



d I ik x 5v y — ik y Sv x \ c 2 s k y 
~dt \ lT J ~ iZpEj 

The asymptotic time dependence of the perturbed potential vorticity can be obtained by integrating 



19 Notice that this is the same time dependence obtained bv lEliassen et aD il953T> in a modal analysis; see the 
discussion surrounding equation 13.611 . These power law time-dependences can be obtained more efficiently by 
solving the large- f limit of equation 13.661 . 
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equation (|3.81[) : 

<Jf ~ f 2 ^ ~ t 1 -™ (3.82) 

for f > 1 and |Ri| -C 1. As noted in §2, an entropy gradient is not required to generate vorticity. 
For N% = 0, a = 1 and the perturbed potential vorticity grows linearly with time. The unstratified 
shearing sheet is recovered in the limit of zero stratification (1/Lp — > 0), since in this limit equation 
(j3.81|) reduces to £ = constant. 

3.5.4 Energetics of the Incompressive Shwaves 

For a physical interpretation of the incompressive shwaves in the stratified shearing sheet, we 
repeat the analysis of section 3.3 for the solution given in the previous section. For a complete 
description of the energy in this case, however, we must include the potential energy of a fluid 
element displaced in the radial direction. Following iMilesl ([196 11 ). an expression for the energy 
in the Boussinesq approximation is obtained by summing equation (|3.46|) multiplied by 5v x and 
equation Q3.47[) multiplied by 5v y . Replacing <5£/£ by £ X /Ls via equation (|3.69f) results in the 
following expression for the energy evolution: 

'' fl V Q 5v 2 + )--E Q Nl£ \ =V Sv x 5v v , (3.83) 



df df \2 2 

where 5v 2 = 5v x + 5vy. The three terms in equation (|3.83j) can be identified as the kinetic energy, 
potential energy and Reynolds stress associated with an individual shwave. One may readily verify 
that the vortical shwaves (see equations (|3.24l) - (|3.26|0 in the unstratified shearing sheet (iVj = 0) 
satisfy equation ()3.83|) . 

The right hand side of equation ()3.83|) can be rewritten — f8v x and individual trailing shwaves 
(f > 0) are therefore associated with a negative angular momentum flux. If the energy were positive 
definite this would require that individual shwaves always decay. But when N% < (Ri < 0) the 
potential energy associated with a displacement is negative, so the energy E\~ can be negative and 
a negative angular momentum flux is not enough to halt shwave growth. 

Our next step is to write the constants of integration C\ and C2 in terms of the initial radial 
velocity and displacement of the shearing wave, 5v x q and £ x o : 

c ^ = qttRi5v x2 {f ) £ x0 + &Ei(fo) 5v x0 ^ = -qQRi 5v x i (f ) g x0 + Ri ^2(^0) Sv x0 ^ 
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where To = k x o/k y , 5v x i is the hypergeometric function given by equation (|3.66f) with C\ = 1 
and C2 = 0, and the other functions are similarly defined. These expressions can be simplified by 
noticing that the denominator of C\ and C 2 is the Wronskian of the differential equation for £, x '- 20 



(l + f 2 )^ + 2f^ + Ri£ E = 0. 

d,T z (XT 



(3.85) 



The Wronskian of this equation is 



w _ d £x2 p 



df 



£ x2 = exp 



2t 2 



1 + T 2 



dr 



l + f a 



(3.86) 



We further simplify the analysis by setting the initial displacement £ x o to zero. 

With these simplifications, the solution given by equations ()3.66|) and (|3.67|) becomes 



5vi 



5v x 



xO 



3 — a 3 + a 1 



1 — a 1 + a 1 ~2i L - 
"T~ ; 2 ; " T ° I 4 ' 4 '2 



, 3 — a 3 + a 3 9 \ /" 5 — a 5 + a 3 > 

Rif F( -^,^^;-;-f 2 fF( — ;-;-f 2 



4 '2 



4 ' 4 '2' 



(3.87) 



5 c 



,r0 



;i+-o 2 ) 



7b f 



3 — a 3 + a 3 



4 '2' T ° 



fn Z F 



1 — a 1 + a 1 



~2 



4 - ; 2 ; - T ) + 



1 — a 1 + a 1 



— ,2' '0 



fpf fF 



3 — a 3 + a 3 



4 '2' 



(3.88) 



As in section 3.3, the energy integral for the incompressive perturbations is given by 



1. 



(E-) = -s F^y k dk (5vi(k )) j desire 



2tt 



Sv x0 



(3.89) 



for initial perturbations perpendicular to and isotropic in k$. Changing integration variables to 
f = q£lt + cot 6, the angular integral becomes 



f 00 
2 / df 

-00 



1 + H - W<5^i(f) + Ri&^f - qnt)5v X 2(f)Y + 



20 Based upon the relationship between a hypergeometric function and its derivatives, 5v x \ 
RiSv x2 = -d(£ x i)/df. 



(3.90) 

d(£ X 2)/dr and 



OS 



where we have used the relation sin# = (1 + Tq) _1 . In the limit of large qflt, the dominant 
contribution to the angular integral comes from the region < f < qVlt. This can be seen from 
the following argument. Using the expansions given by equations (j3.73|) and (|3.76j) . we find the 
angular integrand is 



2\f(f - qVLt)^' 1 [(YiXx + sgn(f)sgn(f - qQt)mV 2 X 2 ) 2 + RLY 2 X 2 (sgn(f ) - sgn(f - q^lt))' 



(3.91) 



for \f \ 1 and |f — q£lt\ ^> 1. Using the relation T(n + 1) = nT(n), one can easily show that 



X 2 = ~^—Vi and V 2 = : -Xi. 



a — 1 



a + 1 



(3.92) 



The integrand therefore simplifies to 



|f(f - qnt^VfXf- [sgn(f) - sgn(f - qUt)f 

1 — a 



(3.93) 



which is zero unless < f < qflt (for t > 0). For large q£lt, therefore, the angular integral is 
approximately given by 



l - a ./„ 



16V? Xf (fqnt) a / f 

— - — — — F a 1 — a- 1 + ar 

a(l - a) gflt V 



qQ,t—v 



, (3.94) 



where 1 <C ^ <C <^7i. For gilt ^ z/, the above expression can be approximated by evaluating it at 
f = qfii, giving 



(Eiiqflt > 1)) ~ leVfXf 



2v2 T(l + a)T(a) ^ = 



a(l-a)r(2a) 



(3.95) 



where we have used equation 15.3.7 in 



Abramowitz and Stegunl (|l972l ) to evaluate -F(a, 6; c; l). 21 
Notice that there is no power-law growth in the perturbation energy for Ri > 1 /4, 22 consistent 
with the classical Richardson criterion (|3.6U|) . In our analysis the energy decays with time for 
2a — 1 < 0, or Ri > 3/16. Thus the energy of an initial isotropic set of incompressive perturbations 
in a radially-stratified shearing sheet-model grows asymptotically (for Ri < 3/16), just like the 

compressive shwaves and unlike the incompressive shwaves in an unstratified shearing sheet, for 

21 We have numerically integrated the angular integral (13.9UI and found this to be an excellent approximation at 
late times. 



For Ri > 1/4, a is imaginary and Ke[t 



t' 1 cos(2|a| In t). 
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which the energy is constant in time. 

The growth of an ensemble of incompressive shwaves in a stratified disk is not due to a Rayleigh- 
Taylor or convective type instability. There is asymptotic growth for < Ri < 3/16, and convective 
instability requires Ri < 0. One can also see this by examining the asymptotic energy for small 
values of |Ri|, such as would be expected for a Keplerian disk with modest radial gradients: 

(Ei(qnt > 1)) ~ [2yr 2 Ri + 0(Ri 2 )] qnt 1 ' 4 ^ ^ (E t (t = 0)). (3.96) 

Evidently for small values of Ri the near-linear growth in time of the energy is independent of the 
sign of Ri and therefore -/V 2 . 23 

3.6 Implications 

We have studied the nonaxisymmetric linear theory of a thin, radially-stratified disk. Our findings 
are: (i) incompressive, short-wavelength perturbations in the unstratified shearing sheet exhibit 
transient growth and asymptotic decay, but the energy of an ensemble of such shwaves is constant 
with time; (ii) short-wavelength compressive shwaves grow asymptotically in the unstratified shear- 
ing sheet, as does the energy of an ensemble of such shwaves, which in the absence of any other 
dissipative effects (e.g., radiative damping) will result in a compressive shwave steepening into a 
train of weak shocks; (iii) incompressive shwaves in the stratified shearing sheet have density and 
azimuthal velocity perturbations <5£, 5v y ~ t~ Rl (for |Ri| <C 1); (iv) incompressive shwaves in the 
stratified shearing sheet are associated with an angular momentum flux proportional to —k x /k y ; 
leading shwaves therefore have positive angular momentum flux and trailing shwaves have negative 
angular momentum flux 24 ; (v) the energy of an ensemble of incompressive shwaves in the stratified 
shearing sheet behaves asymptotically as t 1_4Rl for |Ri| <C 1. For Keplerian disks with modest 
radial gradients, |Ri| is expected to be <C 1, and there will therefore be weak growth in a single 
shwave for Ri < and near-linear growth in the energy of an ensemble of shwaves, independent of 
the sign of Ri. 

Along the way we have found the following solutions: (i) an exact solution for nonvortical 
shwaves in the unstratified shearing sheet, equations (|3.14f) . (|3.15|) and 1)3.20(1 : (ii) a WKB-solution 

23 This asymptotic expression assumes Ri 7^ 0. Notice that the energy at late times can have the opposite sign to 
the initial energy because the potential energy is negative for N x < 0. 

This is consistent with the asymptotic result one obtains from a WKB analysis of incompressive waves fealbusl 

I2003D . 
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for the nonvortical, compressive shwaves in the short- wavelength, high-frequency limit, equations 
(|3.34|) - (|3,36j) : (iii) a solution for incompressive shwaves in the unstratified shearing sheet valid in 
the short- wavelength, low-frequency limit, equations (|3.24j) - ([3.26j) : (iv) a solution for incompressive 
shwaves in the radially-stratified shearing sheet (also valid in the short-wavelength, low-frequency 
limit), equations (ISTofill - (137701) . 

Our results are summarized in Figure l3~T( which shows the regions of amplification and decay 
for shwaves in a stratified disk in the N^/fl, q plane. 

The presence of power-law growth of incompressive shwaves in stratified disks opens the pos- 
sibility of a transition to turbulence as amplified shwaves enter the nonlinear regime. Any such 
transition would depend, however, on the nonlinear behavior of the disk after the shwaves break. 
It is far from clear that they would continue to grow. We will evaluate the nonlinear behavior of 
the disk in subsequent work. 

Our results are essentially in agreement with the numerical results presented bv iKlahrl ((2004), 



that is, we find that arbitrarily large amplification factors can be obtained by starting with appro- 
priate initial conditions. Our results, however, clarify the nature and asymptotic time dependence 
of the growth. O ur r esults on the un s tratifi ed shearing sheet are also consistent with the results of 



rd 



Shepherd! (|1985T ) and 



Afshordi et al 



(2004]), who find that an isotropic ensemble of incompressive 



shwaves have fixed energy. 
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Figure 3.1 A summary of analytic results for shwaves (shearing waves) in a stratified disk. The 
relevant parameters are the local dimensionless shear rate q = — ^dlnfl 2 /dlnr and the dimension- 
less Brunt- Vaisala frequency N 2 /Q 2 . The expected location of a thin, smooth disk is shown as a 
vertically extended ellipse near q = 1.5, N 2 /fl 2 = 0. The far right region (shaded in the figure) is 
forbidden by the H0iland criterion. When q = shear is absent and a modal analysis is possible; 
instability is present for N 2 < 0. Solitary shwaves with Ri = N 2 /(q 2 ft 2 ) < experience asymptotic 
power-law growth (oc for small Ri); since each shwave grows the energy of an ensemble of 
shwaves does as well. For < Ri < 3/16 solitary shwaves decay but the energy of an ensemble of 
shwaves grows as a power-law in time. For Ri > 3/16 both solitary shwaves and the energy of an 
ensemble of shwaves asymptotically decay. 
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Nonlinear Stability of Thin, 
Radially-Stratified Disks 



4.1 Chapter Overview 



We perform local numerical experiments to investigate the nonlinear stability of thin, radially- 
stratified disks. We demonstrate the presence of radial convective instability when the disk is 
nearly in uniform rotation, and show that the net angular momentum transport is slightly inwards, 
consistent with previous investigations of vertical convection. We then show that a convectively- 
unstable equilibrium is stabilized by differential rotation. Convective instability (corresponding to 
Ri — > — oo, where Ri is the radial Richardson number) is suppressed when Ri > —1, i.e. when the 
shear rate becomes greater than the growth rate. Disks with a nearly-Keplerian rotation profile 
and radial gradients on the order of the disk radius have Ri > —0.01 and are therefore stable to 
local nonaxisymmetric disturbances. One implication of our results is that the "Global Baroclinic 



Instability" claimed by 



Klahr and Bodenheimerl (|2003h is either global or nonexistent. 



4.2 Introduction 



In order for astrophysical disks to accrete, angular momentum must be removed from the disk ma- 
terial and transported outwards. In many disks, this outward angular momentum transport is likely 
mediated internally by magnetohydrodvnamic (MHD) turbulence driven by the magnetorotational 



instability (MRI; see 



Balbus and HawlevBl998f ) . A key feature of this transport mechanism is that it 



arises from a local shear instability and is therefore very robust. In addition, MHD turbulence trans- 
ports angular momentum outwards; some other f orms of turbulence, such as convective turbulence, 



appear to transport angular momentum inwards (St one and Balbusl . 11996). The mechanism is only 
effective, however, if the plasma in the disk is sufficiently ionized to be well-coupled to the magnetic 
field (see 31.3.2j) . In portions of disks around young , low-mass stars, in cataclysmic- variable disks 



in quiescence , and in X-ray transients in quiescence ((Stone et al 



Menou . 



2000; 



Gammie and Menou 



1998 



2000), the plasma may be too neutral for the MRI to operate. This presents some difficul- 
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ties for understanding the evolution of these systems, since no robust transport mechanism akin to 



MRI-induced turbulence has been established for p urely-hydrodynamic Kep 



erian shear flows. 



Such a mechanism has been claimed recently bv lKlahr and Bodenheimerl (|2003l ) . who find vor- 
tices and an outward transport of angular momentum in the nonlinear outcome of their global 
simulations. The claim is that this nonlinear outcome is due to a local instability (the "Global 
Baroclinic Instability" ) resulting from the pres ence of an eq uilibrium entropy gradient in the radial 



direction. The instability mechanism invoked ( Klahrl. 12004 ) is an interplay between transient am- 
plification of linear disturbances and nonlinear effects. The existence of such a mechanism would 
have profound implications for understanding the evolution of weakly-ionized disks. 

In Chapter El we have performed a linear stability analysis for local nonaxisymmetric distur- 
bances in the shearing- wave formalism. While the linear theory uncovers no exponentially- growing 
instability (except for convective instability in the absence of shear), interpretation of the results 
is somewhat difficult due to the nonnormal nature of the linear differential operators 1 : one has a 
coupled set of differential equations in time rather than a dispersion relation, which results in a 
nontrivial time dependence for the perturbation amplitudes S(t). In addition, transient amplifica- 
tion does occur for a subset of initial perturbations, and linear theory cannot tell us what effect this 
will have on the nonlinear outcome. For these reasons, and in order to test for the presence of local 
nonlinear instabilities, we here supplement our linear analysis with local numerical experiments. 

We begin in §2 by outlining the basic equations for a local model of a thin disk. In §3 we 
summarize the linear theory results from Chapter|3J We describe our numerical model and nonlinear 
results in §§4 and 5, and discuss the implications of our findings in §6. 



4.3 Basic Equations 



The simulations of 



Klahr and Bodenheimerl ([20031 ) are two-dimensional (without vertical structure) , 



since the salient feature supposedly giving rise to the instability is a radial entropy gradient. The 
simplest model to use for a local verification o f their global results is the two-dimensional shearing 



sheet (see, e.g., 



Goldreich and Tremainc 



1978). This local approximation is made by expanding the 



equations of motion in the ratio of the disk scale height H to the local radius R, and is therefore only 

valid for thin disks (H/R -C 1). The vertical structure is removed by using vertically-integrated 
A nonnormal operator is one that is not self-adjoint, i.e. it does not have orthogonal eigenfunctions. 
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quantities for the fluid variables 2 . The basic equations that one obtains are 



dv VP 

— + — — + 2fl x v - 2qn 2 xx = 0, (4.2) 

Ctt Z-i 



d InS , t , 

"ST = °- <°> 

where £ and P are the two-dimensional density and pressure, S = is the fluid entropy, 3 v 

is the fluid velocity and d/dt is the Lagrangian derivative. The third and fourth terms in equation 
(|4.2j) represent the Coriolis and centrifugal forces in the local model expansion, where Q is the 
local rotation frequency, x is the radial Cartesian coordinate and q is the shear parameter (equal to 
1.5 for a disk with a Keplerian rotation profile). The gravitational potential of the central object 
is included as part of the centrifugal force term in the local-model expansion, and we ignore the 
self-gravity of the disk. 

4.4 Summary of Linear Theory Results 

An equilibrium solution to equations (j4,lj) through l|4 .(-{)) is 

P = P (x), (4.4) 



E = So(ar), (4.5) 



v = v = (-qUx+ ij^) V, ( 4 -6) 

2 This vertical integration is not rigorous; we are assuming that important vertical structure does not develop to 
affect our results. 

3 F or a non-self-gravitating disk the two-dimensional adiabatic index 7 = (37313 — 1)/ (730 + 1) (e.g. lOoldreich et alJ 
Il98fj) . 
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where a prime denotes an x derivative. One can regard the background flow as providing an effective 
shear rate 

qn = -v' (4.7) 

that varies with x, in which case vq = — J x q(s)ds Qy. Due to this background shear, localized 
disturbances can be decomposed in terms of "shwaves", Fourier modes in a frame comoving with 
the shear. These have a time-dependent radial wavenumber given by 

k x (t, x) = k x0 + q(x)Q.k y t. (4.8) 

where k X Q and k y are constants. Here k y is the azimuthal wave number of the shwave. 
In the limit of zero stratification, 

Po(x) — > constant, (4-9) 

T,q(x) — > constant, (4-10) 

v -> -qtlxy, (4.11) 

q^q, (4.12) 

and 

kx -> k x = k x0 + qttkyt. (4.13) 

In Chapter |5J we analyze the time dependence of the shwave amplitudes for both an unstratified 
equilibrium and a radially-stratified equilibrium. As discussed in more detail in Chapter|31 applying 
the shwave formalism to a radially-stratified shearing sheet effectively uses a short-wavelength WKB 
approximation, and is therefore only valid in the limit k y L S> 1, where the background varies on a 
scale L ~ H <C R. The disk scale height H = c s £l, where c s = \/^Pq/Y1q. 

There are three nontrivial shwave solutions in the unstratified shearing sheet, two nonvortical 
and one vortical. The radial stratification gives rise to an additional vortical shwave. In the limit 
of tightly- wound shwaves {\k x \ ^> k y ), the nonvortical and vortical shwaves are compressive and 
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incompressive, respectively. The former are the extension of acoustic modes to nonaxisymmetry, 
and to leading order in (k y L)~ l they are the same both with and without stratification. Since 
the focus of our investigation is on convective instability and the generation of vorticity, we repeat 
here only the solutions for the incompressive vortical shwaves and refer the reader to Chapter |3] for 
further details on the nonvortical shwaves. 

In the unstratified shearing sheet, the solution for the incompressive shwave is given by: 

Sv x = 5v x0 ^, (4.14) 



Sv,, 



kx 



5v x 



(4.15) 



and 



6P 
7P0 



ic s ky 



£*!*+2( g -l)n** 

Ky Cg C S 



(4.16) 



where k 2 = k x + ky, (ko,Sv x o) are the values of (k,5v x ) at t = and an overdot denotes a time 
derivative. 4 

The kinetic energy for a single incompressive shwave can be defined as 



1 



E k = -E (Sv x + 5v. 



(4.17) 



an expression which varies with time and peaks at k x = 0. If one defines an amplification factor 
for an individual shwave, 



A 



E k (k x = 0) 



k 2 

E k (t = 0) + J^> 



(4.18) 



it is apparent that an arbitrary amount of transient amplification in the kinetic energy of an 
individual shwave can be obtained as one increases the amount of swing for a leading shwave 
(k x0 < -k y ). 

This transient amplification of local nonaxisymmetric disturbances is reminiscent of the "swing 
amplification" mec hanism that occurs in disks that are marginally-stable to the axisymmetric gravi- 



tatio nal instability ( Goldreich and Lvnden-Belll 



1965 



Julian and Toomre. 



1966; 



Goldrei ch and Trema inc 



1978). In that context, nonaxisymmetric shwaves experience a short period of exponential growth 



4 As discussed in Chapter |3J this solution is valid for all time only in the short- wavelength limit (k y H 3> 1); for 
Hk y < O(l), an initially-leading incompressive shwave will turn into a compressive shwave near k x = 0. 
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near k x = as they swing from leading to trailing. In order for this mechanism to be effective 
in destabilizing a disk, however, a feedback m echanism is required to convert trailing shwaves into 



leading shwaves ([Binnev and Tremaine 



1987). The arbitrarily-large amplification implied by equa- 



tion (|4.18|) has led so me authors to argue for a bypass transition to turbu l ence in hydrodynami c 



Keplerian shear flows (|Chagelishvili et al 



2003; 



Umurhan and Regev, 



2004; 



Afshordi et al 



2004) 



The reasoning is that nonlinear effects somehow provide the necessary feedback. We show in Chap- 
ter 01 that a ensemble of incompressive shwaves drawn from an isotropic, Gaussian random field 
has a kinetic energy that is a constant, independent of time. This indicates that a random set of 
vortical perturbations will not extract energy from the mean shear. It is clear, however, that the 
validity of this mechanism as a transition to turbulence can only be fully explored via numerical 
experiments. No numerical experiments to date have demonstrated a transition to turbulence in 
Keplerian shear flows. 

In the presence of radial stratification, there are two linearly-independent incompressive shwaves. 
The radial- velocity perturbation satisfies the following equation (we use a subscript s to distinguish 
the stratified from the unstratified case): 



„ 2 d 2 5v xs d5v xs 

(1 + r j — —75 h 4r — — h (Ri + 2)6i 

dr A dr 



>v xs = 0, 



where 



is the time variable and 



f = k x /k y = qttt + k xQ /ky 



Ri 



Nl 



my 



(4.19) 



(4.20) 



(4.21) 



is th e Richardson number, a measur e of the relative importance of buoyancy and shear ((Mile? 



1961 



Howar 



1961; 



Chimona; 



197ft ) 5 . Here 



Ni 



L S L P 



(4.22) 



is the square of the Brunt- Vaisala frequency in the radial direction, where Lp = jPq/Pq and 

Ls = 75*0/5*0 are the equilibrium pressure and entropy length scales in the radial direction. The 

5 As discussed in Chapter [3] equation 14.191 is the same equation that one obtains for the incompressive shwaves 
in a shearing, stratified atmosphere. 



78 



solutions for the other perturbation variables are related to 5v xs by 



Svy S = —f6v xs , (4.23) 



^~ I 5v xs dt (4.24) 



and 

5P S [_„ d ( Sv xs \ : n ,_ ^Sv xs ] {425) 



Pq ic s ky 



d ( 5v xs \ , , 5v x 



df 



Since the solutions to equation Q4.19JI are hyper geometric functions, which have a power-law time 
dependence, it cannot in general be accurately treated with a WKB analysis; there is no asymptotic 
region in time where equation (|4.19|) can be reduced to a dispersion relation. If, however, there 
is a region of the disk where the effective shear is zero, f — ► constant and equation (|4.19f) can be 
expressed as a WKB dispersion relation: 

K x0 S 

with 5(t) oc exp(— iujt). For q ~ and iVj < 0, then, there is convective instability. For disks with 
nearly-Keplerian rotation profiles and modest radial gradients, q ~ 1.5 and one would expect that 
the instability is suppressed by the strong shear. Due to the lack of a dispersion relation, however, 
there is no clear cutoff between exponential and oscillatory time dependence, and establishing a 
rigorous analytic stability criterion is difficult. 

For q y^z 0, the asymptotic time dependence for each perturbation variable at late times is 

6P 8 ex 6v x8 ~ t*** , (4.27) 



<5S S oc 6v ys ~ i 2 ^, (4.28) 

where 



= Vl -4Ri. (4.29) 



The density and y- velocity perturbations therefore grow asymptotically for a > 1, i.e. Ri < 0. For 
small Richardson number, as is expected for a nearly-Keplerian disk with modest radial gradients, 
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a ~ 1 — 2Ri and the asymptotic growth is extremely slow: 

5Z S ~ 5v ys ~ t" Ri . (4.30) 

The energy of an ensemble of shwaves grows asymptotically as t 2a ~ l , or t 1_4Rl for small Ri. The 
ensemble energy growth is thus nearly linear in time for small Ri, independent of the sign of Ri. 6 

The velocity perturbations are changed very little by a weak radial gradient. One would there- 
fore expect that, at least in the linear regime, transient amplification of the kinetic energy for an 
individual shwave is relatively unaffected by the presence of stratification. There is, however, an 
associated density perturbation in the stratified shearing sheet that is not present in the unstratifed 
sheet. 7 This results in transient amplification of the potential energy of an individual shwave. We 
do not derive in Chapter [21 1 a general closed-form expression for the energy of an ensemble of 
incompressive shwaves in the stratified shearing sheet, so it is not entirely clear what effect this 
qualitatively new piece of the energy will have on an ensemble of shwaves in the linear regime. 

In any case, the question of whether or not radial stratification can play a role in generating 
turbulence by interacting with the transient amplification of linear disturbances or by some other 
nonlinear mechanism can only be fully answered with a nonlinear study. For this reason, and due 
to the subtleties involved in the linear analysis, we now turn to the main focus of this paper, which 
is a series of local numerical experiments in a radially-stratified shearing sheet. 



4.5 Numerical Model 



To investigate local nonlinear effects in a radially-stratified thin disk, we int egrate the governing 
equat ions (|4.1|) through Q4.3JI with a hydrodynamics code based on ZEUS ((Stone and Norman . 
1992). This is a time-explicit, operator-split, finite-difference method on a staggered mesh. It uses 
an artificial viscosity to capture shocks. The computational grid is L x x L y in physical size with 
N x x N y grid cells, where x is the radial coordinate and y is the azimuthal coordinate. The boundary 
conditions are periodic in the y-direction and s hearing-period i c in th e x-direction. T he shearing-box 



bou ndary condition s are described in detail in 



and 



Ha wlev et al 



(1223). A s described in lMassetl |2000) 



Gammid (|200lh . advection by the linear shear flow can be done by interpolation. Rather than 



6 We show in Chapter |3] that there is also linear growth in the energy of an ensemble of compressive shwaves. 
7 The amplitude of the density perturbation in the unstratified sheet is an order-of-magnitude lower than the 
velocity perturbations in the short-wavelength limit; see equation l4.16t . 
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using a linear interpolation scheme as in Gammie (2001), we now do the shear transport with the 
same upwind advection algorithm used in the rest of the code. This is less diffusive than linear 
interpolation, and the separation of the shear from the bulk fluid velocity means that one is not 
Courant-limited by large shear velocities at the edges of the computational domain. 

We use the following equilibrium profile, which in general gives rise to an entropy that varies 
with radius: 



ho(x) = h a 



1 — e cos 



2irx 



E (x) 



Mr - 1) 

TK 



i 

r-i 



, P (x)=K^, 



(4.31) 



where h a , e, T and K are model parameters. The flow is maintained in equilibrium by setting the 
initial velocity according to equation (|4.6|) . This equilibrium yields a Brunt- Vaisala frequency 



,2 



N 2 x (x) 



(7 -IX 

7 (r-i)V 



(4.32) 



which can be made positive, negative or zero by varying 7 — V. 

We fix some of the model parameters to yield an equilibrium profile that is appropriate for a 
thin disk. In particular, we want H/Lp ~ H/R <C 1 in order to be consistent with our use of 
a razor-thin (two-dimensional) disk model. In addition, we want the equilibrium values for each 
fluid variable to be of the same order to ensure the applicability of our linear analysis. These 
requirements can be met by choosing K = 1, e = 0.1, L x = 12 and h a = c^T/(T — 1), where 



c s = y (Po/Eo) = 1 is (to within a factor of yfy) the x average of the sound speed. Since the 
equilibrium profile changes with T, we choose a fixed value of T = 4/3, which for 7 = T corresponds 
to a three-dimensional adiabatic index of 7/5. These numbers yield \H/Lp\ < 0.2. Our unfixed 
model parameters are thus L y , q and 7. 

The sinusoidal equilibrium profile we are using generates radial oscillations in the shearing sheet 
due to truncation error. We apply an exponential-damping term to the governing equations in order 
to reduce the spurious oscillations and therefore get cleaner growth-rate measurements. We damp 
the oscillations until their amplitude is equal to that of machine-level noise, and subsequently apply 
low-level random perturbations to trigger any instabilities that may be present. 

As a test for our code, we evolve a particular solution for the incompressive shwaves in the 
radially-stratified shearing sheet (equations |4.19j and |4.23j - |4~23] ). The initial conditions are 



5v x /c s = 5S/Sq = 1 x 10 4 and k. 



'xQ 



-128tt/L x . We set L y = 0.375 and k y = 2n/L y 



m 
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Table 4.1. Summary of Code Runs 



Run 


Description 


Figure(s) 


1 


Linear theory test 


IOIO 


2 


External potential, $7 = 


Pll4"Hl 


q 


External potential, 17 = 1 


14. /'H4.8I 


4 


Uniform rotation (q = 0) 


I4.9B4.11I 


5 


External potential, $7 = 1, boost 


wm 


6 


Small shear (— oo < Ri < —1) 


|4"T3| 


7 


Small shear (q = 0.2, Ri > -1) 


|4.14| 


8 


Aliasing (q = 0.2, Ri > -1) 


14.151 


9 


Parameter survey 




10 


Keplerian disk {q = 1.5, Ri ~ -0.004) 


14.1 VI 



order to operate in the short-wavelength regime, and the other model parameters are q = 1.5 and 
7 — r = —0.3102. The latter value yields a minimum value for N^(x) of —0.01. The results of the 
linear theory test are shown in Figures 14.11 and 14.21 



4.6 Nonlinear Results 

Table H~T1 gives a summary of the runs that we have performed. A detailed description of the setup 
and results for each is given in the following subsections. Our primary diagnostic is a measurement 
of growth rates, and the probe that we use for these measurements is an average over azimuth of 
the absolute value of v x = 5v x at the minimum in N^. Measuring v x allows us to demonstrate 
the damping of the initial radial oscillations, and the average over azimuth masks the interactions 
between multiple WKB modes with different growth rates in our measurements. We will reference 
this probe with the following definition: 

Vt = {\v x (x m in,y)\), (4.33) 

where here angle brackets denote an average over y and x m i n is the x-value at which N^ix) is a 
minimum. 
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Figure 4.1 Evolution of the radial velocity amplitude for a vortical shwave in the radially-stratified 
shearing sheet (Run 1). The heavy line is the analytic result, and the light lines are runs with 
a numerical resolution of (in order of increasing accuracy) N x x N y = 1024 x 16, 2048 x 32 and 
4096 x 64. The number of grid cells are chosen so that the shwave initially has the same number 
of grid cells per wavelength in both the x and y directions. The results are shown for a test point 
at the minimum in N%. 
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Figure 4.2 Evolution of the density amplitude for a vortical shwave in the radially-stratified shearing 
sheet (Run 1). The heavy line is the analytic result, and the light lines are runs with a numerical 
resolution of (in order of increasing accuracy) N x x N y = 1024 x 16, 2048 x 32 and 4096 x 64. 
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4.6.1 External Potential in Non- Rotating Frame 

As a starting problem, we investigate a stratified flow with vq = 0. Such a flow can be maintained 
in equilibrium by replacing the tidal force in equation Q4.2J1 with an external potential <I> = — ho. 
This can be done in either a rotating or non-rotating frame. It is a particularly simple way of 
validating our study of convective instability in the shearing sheet. The condition vq = implies 
q = 0, and therefore equation (|4.26|) should apply in the WKB limit, with the expected growth 
rate obtained by evaluating equation (|4.32fl locally. 8 We have performed a fiducial run with an 
imposed external potential in a non-rotating frame (0 = in equation |4.2j ) to compare with the 
outcome expected from the Schwarzschild stability criterion implied by equation (|4.26|) . We set 
7— T = —0.3102, corresponding to N% min = —0.01, and L y = L x . The expected growth rate for this 
Schwarzschild-unstable equilibrium is 0.1 (in units of the average radial sound-crossing time). The 
numerical resolution for the fiducial run is 512 x 512, and all the variables are randomly perturbed 
at an amplitude of 1.0 x 10~ 12 . 

A plot of vt as a function of time is given in Figure I4.3[ showing the initial damping followed 
by exponential growth in the linear regime. The analytic growth rate is shown on the plot for 
comparison. A least-squares fit of the data in the range 100 < t < 250 yields a measured growth 
rate of 0.0978. 9 Figure B~4"l shows a cross section of N% as a function of x after the instability has 
begun to set in, and Figure l4~5l shows cross sections of the entropy early and late in the nonlinear 
regime. The growth is initially concentrated near the minimum points in N^{x). Eventually 
the entropy turns over completely and settles to a nearly constant value. Figure 14.61 shows two- 
dimensional snapshots of the entropy in the nonlinear regime. Runs with the same equilibrium 
profile except 7 — V > are stable. There is also a long- wavelength axisymmetric instability that is 
present for 7 — T < even in the absence of the small-scale nonaxisymmetric modes. We measure 
its growth rate to be 0.07. Due to the long-wavelength nature of these modes, they are not treatable 
by a local linear analysis. 

4.6.2 External Potential in Rotating Frame 

We have performed the same test as described in §5.1 in a rotating frame (O = 1 in equation 

|4.2j ). Figure l4~7l shows the exponential growth in the linear regime for this run, with a measured 

8 The fastest growing WKB modes will be the ones with a growth rate corresponding to the minimum in N%. 
9 Measurements of the growth rate earlier in the linear regime or over a larger range of data yield results that differ 
from this value by at most 5%. 
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growth rate of 0.0977. Figure 14.81 shows snapshots of the entropy in the nonlinear regime. The 
results are similar to the nonrotating case, except that 1) rotation suppresses the long-wavelength 
axisymmetric instability; 2) the nonlinear outcome exhibits more coherent structures in the rotating 
case including transient vortices; and 3) these coherent structures eventually become unstable to a 
Kelvin-Helmholz-ty p e inst ability. 

4.6.3 Uniform Rotation 

Having demonstrated the viability of simulating convective instability in the local model, we now 
turn to the physically-realistic equilibrium described in §4. We begin by setting the shear parameter 
q to zero in order to make contact with the results of §§5.1 and 5.2. This is analogous to a disk in 
uniform rotation. The other model parameters are the same as for the previous runs. While there 
is still an effective shear —0.05 < q < 0.05, near q = one expects the modes to obey equation 
(j4.26|) in the WKB limit. Figures 14.91 and 14.101 give the linear and nonlinear results for this run. 
The measured growth rate in the linear regime is 0.0809. 



C onsistent wit h results from numerical simulations of vertical convection ((Stone and Ba lbu 



1996 



Cabot 



1996]), the angular momentum transport associated with radial convection is inwards. 



Figure 14.111 shows the evolution of the dimensionless angular momentum flux 

a = t 1 ,tj\ / ^5v x 5v y dxdy, (4.34) 



where {Pq} is the radial average of the equilibrium pressure, for an extended version of Run 4. 
Averaging over the last 1200O" 1 yields a ~ —10 . 

There are two reasons for the larger error in the measured growth rate for this run: 1) the 
equilibrium velocity gives rise to numerical diffusion due to the motion of the fluid variables with 
respect to the grid; and 2) since the growing modes are being advected in the azimuthal direction, 
the maximum growth does not occur at the grid scale. The latter effect can be seen in Figure H. 101 
several grid cells are required for a well-resolved wavelength. In order to resolve smaller wavelengths, 
we have repeated this run with L y = 6, 3 and 1.5. The results are plotted in Figure H31 along with 
the results from the L y = 12 run. The measured growth rate for the L y = 1.5 run is 0.0924. 

To quantify the effects of numerical diffusion, we have performed a series of tests similar to 
Run 2 (external potential in a rotating frame) but with an overall boost in the azimuthal direc- 
tion. Figure H. 121 shows measured growth rates as a function of boost at three different numerical 



resolutions. The largest boost magnitude in this plot corresponds to the velocity at the minimum 
in N x for a run with q = 1.5. 

4.6.4 Shearing Sheet 

To investigate the effect of differential rotation upon the growth of this instability, we have per- 
formed a series of simulations with nonzero q. Intuitively, one expects the instability to be sup- 
pressed when the shear rate is greater than the growth rate, i.e. for Ri > —1. Figure H.13I shows 
growth rates from a series of runs with N% min = —0.01 and small, nonzero values of q at three 
numerical resolutions. This figure clearly demonstrates our main result: convective instability is 
suppressed by differential rotation. The expected growth rate from linear theory (ypV|| at q = 0) 
is shown in Figure H. 131 as a dotted line. If there is a radial position where q(x) = (i.e., Ri = — oo), 
vt at that position looks similar to that of the previous runs (very little deviation from a straight 
line); these measurements are indicated on the plot with solid points. For q > 0.055 there is no 
longer any point where q{x) = 0; in that case vt was measured at the radial average between the 
minimum in N x (x) and the minimum in q(x), since this is where the maximum growth occured. 
The data for these measurements, which are indicated in Figure 14.131 with open points, is not as 
clean as it is for the runs with Ri = — oo (see Figure H. 14(1 . All of the growth rate measurements in 
Figure l4~T3l were obtained by a least-squares fit of the data in the range 1 x 10~ 9 < v t /c s < 1 x 10~ 5 . 
The dashed line in Figure R. 131 indicates the value of q for which Ri m j n = — 1. 

Some of the growth in Figure 14.131 appears to be due to aliasing. This is a numerical effect 
in finite-difference codes that results in an artificial transfer of power from trailing shwaves into 
leading shwaves as the shwave is lost at the grid scale. One expects aliasing to occur approximately 
at intervals of 

Af = ^^, (4.35) 

fly l^x 

where n y is the azimuthal shwave number. This interval corresponds to Ak x (t) = 2n/dx, where 
dx = L x /N x is the radial grid scale. Based upon expression (|4,35[) . aliasing effects should be more 
pronounced at lower numerical resolution because the code has less time to evolve a shwave before 
the wavelength of the shwave becomes smaller than the grid scale. It can be seen from the far- 
right data point in Figure H. 131 ( Run 7 in Table l4~Tj) that the measured growth rate decreases with 
increasing resolution. The evolution of x>t for this run is shown in Figure H. 141 

The effects of aliasing can be seen explicitly by evolving a single shwave, as was done for our 
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linear theory test (Figures 14.11 and 14. 2[) . Figure 14.151 shows the evolution of the density pertur- 
bation for a single shwave using the same parameters that were used for Run 7: L x = L y = 12, 
Nxmin = — 0.01 and q = 0.2. The initial shwave vector used was (k X Q, k y ) = (—8TT/L x ,8ir/L y ). This 
corresponds to n y = 4, and the expected aliasing interval ()4.35[) is therefore At = N x /4. Runs at 
three numerical resolutions are plotted in Figure 11.151 and the aliasing interval at each resolution 
is consistent with expression (|4.35fl . It is clear from Figure 11.151 that a lower resolution results in 
a larger overall growth at the end of the run. It also appears that the growth seen in Figure 14.151 
requires a negative entropy gradient. We have performed this same test with N x > 0, and while 
aliasing occurs at the same interval, there is no overall growth in the perturbations. This is likely 
due to the fact that the perturbations decay asymptotically for N x > (see expression |4.30| ). 

Figure H. 161 summarizes the parameter space we have surveyed, indicating that there is insta- 
bility only for q ~ and N x < 0. The numerical resolution in all of these runs is 512 x 512. 
Figure 13 . 1 71 shows the evolution of the radial velocity in Run 10, a run with realistic parameters for 
a disk with a nearly-Keplerian rotation profile and radial gradients on the order of the disk radius: 
q = 1.5 and N^ min = —0.01 (corresponding to Ri ~ —0.004). Clearly no instability is occurring 
on a dynamical timescale. This plot is typical of all runs for which the evolution was stable. To 
give a sense for the minimum growth rate that we are able to measure, we have also plotted in 
Figure II. 17l the results from several unstable runs with q = and a boost equivalent to the velocity 
at the minimum in Nl for Run 10. It is difficult to measure a growth rate for the smallest value of 
^xmini but ^ ^ s clear that there is activity present in this run which does not occur in the stable 
run. Based upon Figure 11.171 a conservative estimate for the minimum growth rate that should be 
detectable in our simulations is 0.00250. 



4.7 Implications 

Our results seem to indicate that nearly-Keplerian disks with weak radial gradients are stable to 
local nonaxisymmetric disturbances, although we cannot exclude instability at very high Reynolds 
number. Figure 14.131 demonstrates that convective instability, present when the shear is nearly 
zero, is stabilized by differential rotation. Perturbations simply do not have time to grow before 
they are pulled apart by the shear. 



An important implication of our results is that the instability claimed by 



Klahr and Bodenheimer 



(2003j) is not a linear or nonlinear local nonaxisymmetric instability. Figure H. 151 suggests that the 
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results of lKlahr and Bodenheime r (2003) may be due, at least in part, to aliasing. They use a finite 
difference code at fairly low numerical resolution (< 128 2 ), and growth is only observed in runs 
with a negative entropy gradient. Curvature effects and the effects of boundary conditions, which 
may also play a role in their global results, cannot be tested in our local model. 




Figure 4.3 Evolution of vt as a function of time for Run 2 (external potential, non-rotating frame). 
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Figure 4.4 Plot of N% (averaged over y) as a function of x for Run 2. The dotted line shows 
the equilibrium profile, and the solid line shows a snapshot during the nonlinear regime. Growth 
initially occurs at the minimum in N%. 




Figure 4.5 Plot of the entropy (averaged over y) as a function of x for Run 2. The dotted line 
shows the equilibrium profile, and the solid lines show snapshots during the nonlinear regime. The 
entropy eventually settles to a nearly-constant value. 
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Figure 4.6 Snapshots of the entropy in the nonlinear regime for Run 2, indicating maximum growth 
for modes near the grid scale and the eventual turnover of the equilibrium entropy profile to its 
average value. Dark shades indicate values above (red in the color version) and below (blue in the 
color version) the average value (yellow in the color version). 
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Figure 4.7 Evolution of vt as a function of time for Run 3 (external potential, rotating frame). The 
dotted line shows the expected growth rate. 
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Figure 4.8 Snapshots of the entropy in the nonlinear regime for Run 3. Dark shades indicate values 
above (red in the color version) and below (blue in the color version) the average value (yellow in 
the color version). 
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Figure 4.9 Evolution of vt as a function of time for Run 4 (q = 0). The dotted line shows the 
expected growth rate, and the solid lines are runs with (in order of increasing growth) L y = 12, 6, 
3 and 1.5. 
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Figure 4.12 Growth rates as a function of azimuthal boost in a series of runs with an external 
potential and N^ min = —0.01. The dotted line shows the analytic growth rate from linear theory. 
The open circle denotes the growth rate that was measured in Run 4, with the boost corresponding 
to the magnitude of the velocity at the minimum in iVj for Run 3 (q = 0). The largest boost 
magnitude corresponds to the velocity at the minimum in iVj for Run 10 (q = 1.5). 
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Figure 4.15 Evolution of the density perturbation for a single shwave with q = 0.2, N^ min = —0.01 
and L y = L x (Run 8). The linear theory result is shown as a dotted line, along with results at three 
numerical resolutions. Aliasing occurs when k x (t) = 2ir/dx. The overall growth, which is greater 
at lower numerical resolution, requires iVj < 0. 
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Figure 4.16 Parameter space surveyed in a search for nonlinear instabilities. Closed (open) circles 
denote runs that were unstable (stable). The only instability found was convective instability for 
q ~ and N% < (Ri — > — oo). (We do not include on this plot the runs shown in Figure H. 131 ) 
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Figure 4.17 Evolution of vt as a function of time for Run 10 (q = 1.5, N^ min = —0.01). Also shown 
are runs with q = and an overall boost equivalent to the velocity at the minimum in Nl for Run 
10, for N% min = —0.01 (measured growth rate of 0.058), N% min = —0.003 (measured growth rate 
of 0.021) and N^ min = -0.001 (measured growth rate of 0.0025). 
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Vortices in Thin, Compressible, 
Unmagnetized Disks 



5.1 Chapter Overview 

We consider the formation and evolution of vortices in a hydrodynamic shearing-sheet model. The 
evolution is done numerically using a version of the ZEUS code. Consistent with earlier results, 
an injected vorticity field evolves into a set of long-lived vortices each of which has radial extent 
comparable to the local scale height. But we also find that the resulting velocity field has positive 
shear stress (E5v r 6vj,) • This effect appears only at high resolution. The transport, which decays 
with time as t~ 1//2 , arises primarily because the vortices drive compressive motions. This result 
suggests a possible mechanism for angular momentum transport in low-ionization disks, with two 
important caveats: a mechanism must be found to inject vorticity into the disk, and the vortices 
must not decay rapidly due to three-dimensional instabilities. 1 



5.2 Introduction 



Astrophysical disks are common because the specific angular momentum of the matter inside them 
is well-conserved. They evolve because angular momentum conservation is weakly compromised, 
either because of diffusion of angular momentum within the disk or because of direct application 
of external torques. 

In astrophysical disks composed of a well-ionized plasma it is likely that some, perhaps most, 
of the evolution is driven by diffusion of angular momentum within the disk. This view is certainly 



consistent with o 



( Baptista et al 



aseryations of steadily accreting cataclysmic variable systems like UX Ursa Majoris 



1998; 



Baptista, 



2004), whose radial surface-brightness profile is consistent with 
steady accretion-flow models in which the bulk of the accretion energy is dissipated within the 
disk. 

Angular momentum diffusion in well-ionized disks is likely driven by magnetohydrodynamic 
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(MHD) turbulence. Analytic analyses, numerical experiments, and laboratory evidence strongly 
suggest that well-couple d plasmas in differentially - rotating; flows are s ubject to the magnetorota- 



tional instability (MRI; 



Balbus and Hawlev 



1991, 



1998; Bal 



But MHD turbulence is 



initia ted by the MRI only so 



field (Kunz and Balbus 



2004; 



long a s the p lasma is sufficiently ionized to couple to the magnetic 



Desch, 



2004). In disks around young stars, cataclysmic- variable and 



neutral to support magnetic activity ( 


Garni) 


trie and Menou. 


1998; 


Menou, 


2000; 


Stone et al. . 


2000: 


Menou and Ouataert 


2001: 


Promang et al.. 


2OO2I1. This motivates interest in non-MHD angular 



momentum transport mechanisms. 

Within the last few years, a body of work has been d eveloped sugge s ting that vortices c an be 
generated as a result of global hydrodynamic ins tability fH awlev 



Haw 



2004 



ev 



1990 



Lovelace et al 



1999 



Li et al 



1987 



Blaes and Hawlev, 



1988 



2000) or local hydrodynamic ins t ability ( Klahr and Bo denheimer 



, that vortices in disks 



Barranco and Marcus 



angular momentum ([Li et al 



may b e long-lived ([Godon and Livid . 1999, 



2000; 



Umurhan and Eeg ev. 



2005,V and that these vortices may be related to an outward flux of 



2001 



Barranco and Marcu; 



2005). If these claims can be verified 



then the consequences for low-ionization disks would be profound. 

Here we investigate the evolution of a disk that is given a large initial vortical velocity perturba- 
tion. Our study is done in the context of a (two-dimensional) shearing-sheet model, which permits 
us to resolve the dynamics to a degree that is not currently possible in a glo bal disk model. Our 



mode 



2004 



is also fully compressible, u nlike previous work using a local model ([Umurhan and Regev . 

Barranco and Marcusl . 2005 ). The former assume incompressible flow and the latter use the 

I I! I 

•ugh 1969), which filters out the high-frequency acoustic waves. We 



anelastic approximation fe.g..lGough 1 
will show that compressibility and acoustic waves play an essential part in the angular momentum 
transport. 

Our paper is organized as follows. In §2 we describe the model. In §3 we describe the evolution 
of a fiducial, high-resolution model. In §4 we investigate the dependence of the results on model 
parameters. And in §5 we describe implications, with an emphasis on key open questions: are 
the vortices destroyed by three-dimensional instabilities?; and do mechanisms exist that can inject 
vorticity into the disk? 
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5.3 Model 



The shearing-sheet model is obtained via a rigorous expansion of the two-dimensional hydrodynamic 
equations of motion to lowest order in H/R, where H = c s /£l is the disk scale height (c s is the 
isothermal sound sp eed and f2 is the local rotation frequency) and R is the local radius. See 



Naravan et al 



( 1987) for a description. Adopting a local Cartesian coordinate system where the x 
axis is oriented parallel to the radius vector and the y axis points forward in azimuth, the equations 
of motion become 

f + EV.. = 0, (5.1) 

dv VP _ ^ 9 A . . 

— + -— + 2fi x v - 2qQ 2 xx = 0, (5.2) 
dt 2-j 

where £ and P are the two-dimensional density and pressure, v is the fluid velocity and d/dt is 
the Lagrangian derivative. The third and fourth terms in equation (|5.2|) represent the Coriolis 
and centrifugal forces in the local model expansion, where q = — (1/2) din Q 2 /dlnr is the shear 
parameter. We will assume throughout that q = 3/2, corresponding to a Keplerian shear profile. 
We close the above equations with an isothermal equation of state 

P = c^S, (5.3) 

where c s is constant in time and space. 

Equations (|5.1jl through (|5.3I) can be combined to show that the vertical component of potential 
vorticity 

€ = (VX " + 2 "^ (5.4) 

is a constant of the motion; i.e., the potential vorticity of fluid elements in two dimensions is 
conserved. 

An equilibrium solution to the equations of motion is 

S = Sq = const. (5-5) 



c^Sq = const. (5-6) 
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(5.7) 



-qQx 



(5i 



Thus the differential rotation of the disk makes an appearance in the f orm of a linear shear . 



We integrate the above equations using a version of the ZEUS code (|Stone and NormanL 11992). 
ZEUS is a time-explicit, operator-split scheme on a staggered mesh. It uses artificial viscosity to 
capture shocks. Our computational domain is a rectangle of size L x x L y containing N x x N y grid 
cells. The numerical resolution is therefore Ax x Ay = L x /N x x L y /N y . 

Our code differs from the standard ZEUS algorithm in two respects. First, we have implemented 
a version of the shearing-box boundary conditions. The model is then periodic in the y direction; 
the x boundaries are initial joined in a periodic fashion, but they are allowed to shear with respect 



to each other, becoming periodic agai n when t = nL „ / (qO L 



n 



of the boundary conditions is given in 



Hawlev et al 



1, 2, A detailed description 



( 199.51 1. 

Second, we treat advection by the mean flow vq = —qVtxy separately from advection by the 
perturbed fl ow 5v _ v — vn- Mean-flow advection can be done by inter polation , using the algorithm 



described in 



Gammid (|200lh . which is similar to the FARGO scheme (Massct 



200(1 ) . This has the 



advantage that the timestep is not limited by the mean flow velocity (it is \5v\ rather than \v\ 
that enters the Courant condition). This permits the use of a timestep that is larger than the 
usual timestep by ~ L x /H if L x S> H. The shear-interpolation scheme also makes the algorithm 
more nearly translation-invariant in the x — y plane, thereby more nearly embodying an important 
symmetry of the underlying equations. 



5.3.1 Initial Conditions 

Without a specific model for the process that is injecting the vorticity, it is difficult to settle on 
a particular set of initial conditions, or to know how these initial conditions ought to vary when 
the size of the box is allowed to vary. Our choice of initial conditions is therefore somewhat 
arbitrary. We use a set of initial (incompressive) velocity perturbations drawn from a Gaussian 
random field. The amplitude of the perturbations is characterized by a = (\5v /cgl 2 ) 1 ^ 2 . The power 
spectrum is \5v\ 2 ~ A; -8 / 3 , corresponding to the energy spectrum (E^ ~ /c~ 5//3 ) of a two-dimensional 
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Kolmogorov inverse turbulent cascade, with cutoffs at k m i n = (1/2){2tt / H) and k max = 32k m i n 2 . 
The surface density is not perturbed. These initial conditions correspond to a set of purely vortical 
perturbations. The parameters for our fiducial run are L x = L y = AH and a = 0.4. 

5.3.2 Code Verification 



Although our basic algorithm has already been tested (see lOammiel l200lh . we test the current 
version of our code by making a comparison with linear theory. Due to the underlying shear, 
small-amplitude perturbations in the shearing sheet are naturally decomposed in terms of shearing 
waves or shwaves (see Chapter |2J), Fourier components in the "co-shearing" frame. These have time- 
dependent wavenumber k(t) = k x (t)x + k y y, where k x {t) = k X Q+qVtk y t and k X Q and k y are constant. 
The evolution of a single Fourier component can be calculated by integrating an ordinary differential 
equation for the amplitude of the shwave. For purely vortical (nonzero potential vorticity) or non- 
vortical perturbations, the evolution can be obtained analytically. The explicit expression for the 
amplitude of a vortical (incompressive) shwave is 

k 2 1 + t 2 

Sv xi = 5v x0 -^ = 5v x0 — — %, (5.9) 

where k 2 = k 2 + fc 2 ,, r = qtlt + k x o/k y and a subscript on a quantity indicates its value at t = 0. 3 
The amplitude of a non-vortical (compressive) shwave satisfies the differential equation 

f^f + (c 2 k 2 + tf) 5v yc = 0, (5.10) 

the solutions of which are parabolic cylinder functions. See Chapter |3] for further details on the 
shwave solutions. 

Figures 15.11 and 15.21 compare the numerical evolution of both vortical and compressive shwave 
amplitudes with their analytic solutions. The initial shwave vector (k x o, k y ) is (— 16ir/L x , 4ir/L y ) for 
the vortical shwave and (—8ir/L x , 2ir/L y ) for the compressive shwave. The other model parameters 
are the same as those in the fiducial run, except that L = 0.5H for the vortical-shwave evolution 
since k y H 3> 1 is required to prevent mixing between vortical and non-vortical shwaves near r = 0. 
The shwaves are well resolved until the radial wavelength = 4 x Ax, and the code is capable of 



2 We have compared our fiducial run to runs with a different range in k, corresponding to vorticity injection either 
at scales ~ H or scales ~ Q.1H . The results are qualitatively the same. 



3 



This solution is valid at all times only for short-wavelength vortical perturbations (kH ^> 1). 
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Figure 5.1 Evolution of the radial velocity amplitude for a vortical shwave. The heavy line is 
the analytic result, and the light lines are numerical results with (in order of increasing accuracy) 
N x = N y = 32, 64, 128 and 256. 



l — i — i — i — | — i — i — i — | — i — i — i — | — i — i — i — p 




Figure 5.2 Evolution of the azimuthal velocity amplitude for a nonvortical shwave. The heavy line 
is the analytic result, and the light lines are numerical results with (in order of increasing accuracy) 
N x = N y = 32, 64, 128 and 256. 
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tracking both potential-vorticity and compressive perturbations with high accuracy. 

Without a specific model for the process that is injecting the vorticity, it is difficult to settle 
on a particular set of initial conditions, or to know how these initial conditions ought to vary 
when the size of the box is allowed to vary. Our choice of initial conditions is therefore somewhat 
arbitrary. We use a set of initial (incompressive) velocity perturbations drawn from a Gaussian 
random field. The amplitude of the perturbations is characterized by a = (|<5i>| 2 ) 1//2 . The power 
spectrum is appropriate for two dimensional Kolmogorov turbulence, \5v\ 2 ~ k~ 8 ^ 3 , with cutoffs 
at k m i n = (l/2)(2ir/H) and k max = 32k m i n . The surface density is not perturbed. These initial 
conditions correspond to a set of purely vortical perturbations. 



5.4 Results 



The evolution of the potential vorticity in our fiducial run is shown in Figure 15.31 The snapshots 
are shown in lexicographic order beginning with the initial conditions, which have equal positive 
and negative 8^. 

One of the most remarkable features of the fiducial run evolution is the appearance of com- 
paratively stable, long-lived vortices. These yortices have nega t ive 6£ and a re theref o re da rk in 
Figur e 15.31 S imilar yortices have been seen by 



Oodon and Livid (1999 



2noovlLi fit, aiJfcnmh and 



Umurhan and Re ecv ( 2004j) . Cross sections of one of the vortices at the end of the run are shown in 
Figures 15.41 and 15.51 In our models the vortices are not associated with easily identifiable features 
in the surface density, since the perturbed vorticity is not large enough to require, through the 
equilibrium condition, an order unity increase in the local pressure. 

While the vortices are long-lived, they do decay. Figure EE1 shows the evolution of the perturbed 
(noncircular) kinetic energy 



1 



E K = -E(5vi + Sv' y ) 



(5.11) 



in the fiducial run. Evidently the kinetic energy decays approximately as t 1//2 (which is remarkable 
in that, if the vortices w ould corresp ond to features in luminosity that decay as t -1 / 2 , they could 



produce nicker noise; see 



Pres; 



1978). Runs with half and twice the resolution decay in the same 
fashion, but if the resolution is reduced to 64 2 the kinetic energy decays exponentially. Resolution 
of at least 128 zones per scale height appears to be required. 

What is even more remarkable is that the vortices are associated with an outward angular 
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Figure 5.3 Panels show the evolution of the potential vorticity in the fiducial run. The size is 
AH x AH and the numerical resolution is 1024 2 . The initial conditions are shown in the upper left 
corner, and the other frames follow in lexicographic order at intervals of 22.2S7 1 . Dark shardes 
(blue and black in the color version) indicate potential vorticity smaller than f2/(2£o); light shades 
(yellow and red in the color version) show positive potential vorticity perturbations. Evidently only 
the "anticyclonic" (negative potential vorticity perturbation) vortices survive. Each vortex sheds 
sound waves, which steepen into trailing shocks. 



110 



Figure 5.4 Radial slice of a vortex at the end of the fiducial run. The heavy line shows the potential 
vorticity, the light line shows the magnitude of the velocity and the dotted line shows the surface 
density. 




Figure 5.5 Azimuthal slice of a vortex at the end of the fiducial run. The heavy line shows the 
potential vorticity, the light line shows the magnitude of the velocity and the dotted line shows the 
surface density. 



Ill 
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momentum flux, due to the driving of compressive motions by the vortices. Figure 15.71 shows the 
evolution of the dimensionless angular momentum flux 

a = t t \ 91 / ^5v x 5v y dxdy (5.12) 
L x LyS c|l J 

for models with a variety of resolutions. The data has been boxcar smoothed over an interval 
At = 10f2 _1 to make the plot readable. Again, a resolution of at least 512 2 appears to be required 
for a converged measurement of the shear stress. For the most highly resolved models a evolves 
like the kinetic energy, oc t -1 / 2 . 




Figure 5.7 Evolution of the shear stress a in the fiducial run and a set of runs at lower resolutions. 

Compressibility is crucial for development of the anglar momentum flux. We have demonstrated 
this in two ways. First, we have taken the fiducial run and decomposed the velocity field into a 
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compressive and an incompressive part (i.e., into potential and solenoidal pieces in Fourier space) 
and measured the stress associated with each. For a set of snapshots taken from the last half of 
the fiducial run, the average total a = 0.0036; the incompressive component is on = —0.0006; the 
compressive component is a c = 0.0032. The remaining alpha a x = 0.00099 is in cro ss-corre l ations 
bet ween the i ncom pressive and compressive pieces of the velocity field. As argued in lBalbusi (|2000) 
and iBalbusl (|2003l ) . both incompressive trailing shwaves and incompressive turbulence tend to 
transport angular momentum inward, whereas trailing compressive disturbances transport angular 
momentum outward. Our negative (positive) value for cti (a c ) is consistent with this. 

Second, we have reduced the size of the model and reduced the amplitude of the initial per- 
turbation so that it scales with the shear velocity at the edge of the model (constant "inten- 
sity" of the turbulence, in Umurhan and Regev's parlance). Thus the Mach number of the tur- 
bulence is reduced in proportion to the size of the box. We have compared four models, with 
L = (4,2, 1,0. 5)H and a = (0.8, 0.4, 0.2, 0.1)c s . We would expect the lower Mach number models 
to have smaller-amplitude compressive velocity fields and therefore, consistent with the above re- 
sults, smaller angular momentum flux a. Averaging over the second half of the simulation, we find 
a = (0.0031,0.0018,7.2 x 10~ 5 ,-9.5 x 10~ 7 ). 

An additional confirmation of our overall picture can be seen in Figure 15.81 in which we show a 
snapshot of the velocity divergence superimposed on the potential vorticity for a medium-resolution 
(256 2 ) version of the fiducial run. 4 The position of the shocks with respect to the vortices in this 
figure is consistent with our interpretation that the former are generated by the latter. 

The smallest of our simulations {L = 0.5H) is nearly incompressible, but we continue to observe 
t~ x l 2 decay (least squares fit pow er law is —0.49) at late times. The reason that we see decay 



while 



Umurhan and Regevl (|200J) do not may be that: (1) the remaining compressibility in our 



model causes added dissipa tion; (2) the numerica" 



Umurhan and Regevl 1)2004 ) : (3) the code used by 



dissipation in our c o de is larger than that of 



Umurhan and Regevl ()2004h could somehow be 



aliasing power from trailing shwaves to leading shwaves (although they do explicitly discuss, and 
dismiss, this possibility). 

To highlight the dangers of aliasing for our finite-difference code, in Figure 15.91 we show the 
evolution of a vortical shwave amplitude at low resolution (64 2 ), in units of r. We use the same 

parameters as those in our linear-theory test (Figures 15.11 and 15.2)1 , for which the initial shwave 

4 At higher resolutions, shocks are generated earlier in the simulation from smaller vortices, and it is more difficult 
to see the effect we are describing due to the random nature of the vortices at this early stage. 
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Figure 5.8 Snapshot of the velocity divergence superimposed on the potential vorticity in a medium- 
resolution (256 2 ) version of the fiducial run. The thin (red in the color version) contours indicate 
negative divergence and are associated with shocks. The thick (blue in the color version) contours 
indicate negative potential vorticity. 
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Figure 5.9 Evolution of a vortical shwave amplitude in a low-resolution (64 2 
The initial shwave vector (k X Q,k y ) is (— 16ir/L x , 4ir/L y ), corresponding to tq 
between successive peaks (a numerical effect due to aliasing) is r = N x /n y , where n 
azimuthal wavenumber. 



run, in units of r. 
= —4. The interval 
2 is the 
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vector corresponds to To = —4. The initially-leading shwave swings into a trailing shwave, the radial 
wavelength is eventually lost near the grid scale, and due to aliasing the code picks up the evolution 
of the shwave again as a leading shwave. Repeating this test at higher resolutions indicates that 
successive swings from leading to trailing occur at an interval of r = N x /n y , where n y = 2 is the 
azimuthal wavenumber of the shwave. This is equivalent to k x (t) = 2ir/Ax. The decay of the 
successive linear solutions with time is due to numerical diffusion. 

Figure 15.91 suggests that it is easier to inject power into the simulation due to aliasing rather 
than to remove power due to numerical diffusion. We do not believe, however, that aliasing is 
affecting our high-resolution results. In addition, if we assume that the flow in our simulations 
can be modeled as two-dimensional Kolmogorov turbulence, then 5v rms ~ A 1 / 3 , where 6v rms is 
the rms velocity variation across a scale A. The velocity due to the mean shear at these scales is 
shear ~ <?OA, and 5v rms /5v shear ~ A~ 2 / 3 . The velocities at the smallest scales are thus dominated 
by turbulence rather than by the mean shear. This conclusion is supported by the convergence of 
our numerical results at high resolution. 

Our model contains two additional numerical parameters: the size L and the initial turbulence 
amplitude a. Figure 15.101 shows the evolution of a for several values of a. Evidently for small 
enough values of a the a amplitude is reduced, but for near-sonic initial Mach numbers the a 
amplitude saturates (or at least the dependence on a is greatly weakened). Figure IB*. 1 II shows the 
evolution for several values of L but the same initial a and the identical initial power spectrum. 
For large enough L the shear stress appears to be independent of L. 

Finally, we have studied the autocorrelation function of the potential vorticity as a means of 
characterizing structure inside the flow. Figures 15.121 and 15.131 shows the autocorrelation function 
measured in the fiducial model and in an otherwise identical model with L = 8H. Evidently the 
potential vorticity is correlated over about one-half a scale height in radius, independent of the size 
of the model. This supports the idea that compressive effects limi t the size of the vortices , since the 



shear flow becomes supersonic across a vortex of size ~ H ((Barge and Sommeria 
200lL 



1995; 



Li et al 



5.5 Conclusion 

The presence of long-lived vortices in weakly-ionized disks may be an integral part of the angular 
momentum transport mechanism in these systems. The key result we have shown here is that 
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Figure 5.10 Evolution of the shear stress a in a set of runs at with varying initial a. Apparently for 
low values of a the shear stress is reduced, but for initial Mach number near 1 the stress saturates. 
All runs have L = AH. 
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Figure 5.11 Evolution of the shear stress a in a set of runs at with varying initial L, but the same 
initial Mach number a. 
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compressibility of the flow is an extremely important factor in providing a significant, positively- 
correlated average shear stress with its associated outward transport of angular momentum. Previ- 



ous results using a 



ocal model have assumed i ncompressible flow and either 



ment um transport (|Umurhan and Regev 



report no angular mo- 



2004) or report a value (a ~ 10 



Barranco and Marcu; 



2005) that is two orders of magnit ude lower than what we find when we inclu de the effects of 



1999, 



2000; 



Li et al 



2001) have a difficult 



compressibility. Global simulations ([Godon and Livic . 
time accessing the high resolution that we have shown is required for a significant shear stress due 
to compressibility. 

Our work leaves open the key question of what happens in three dimensions. Our vortices, 
which have radial and azimuthal extent < H, are in herently three-d imensional. Three-dimensional 
vortices are susceptible to the elliptical instability (|Kerswelll . 120021 ) and are likely to be destroyed 
on a dynamical timescale. The fact that vortices persist in o ur two-dimens i onal s imulations and 
not in the local (three-dimensional) shearing-bo x calculations of 



dimensionality. The recent numerical results of 



Barranco and Marcu; 



Balbus et al. ( 1996) is likely due to 



2005) indicate that vortices 



near the disk midplane are quickly destroyed, whereas vortices survive if they are a couple of scale 
heights away from the midplane. Strong vertical stratification away from the midplane may enforce 



two-dimensional flow and a. 
The initial conditions in 



low the vortices that we consider here to survive. 



Barranco and Marcus (2005) are analytic solutions for two-dimensional 



vortices that are stacked into a three-dimensional column. The stable, off-midplane vortices ap- 
parently arise due to the breaking of internal gravity waves generated by the midplane vortices 
before they become unstable. There is also an unidentified instability that breaks a single off- 
midplane vortex into several vortices. These simulations leave open the question of whether stable 
off-midplane vortices can be generated from a random set of initial vorticity perturbations rather 
than the special vortex solutions that are imposed. 

Our work also leaves open the key question of what generates the initial vorticity. One possibility 
is that material builds up a t part icular radii in the disk, resulting in a global instabili t y (e.g . 



Papaloizou and Pringle 



1984 



1985) and a breakdown of the flow into vortices ([Li et al 



2001). 



Another possibility for vortex generation in variable systems i s that the MHD turbulence , which 



1998ft . leaves 



likely operates during an outburst but decays as the disk cools (jOammie and MenouL 
behind some residual vorticity. The viability of such a mechan ism c ould be tested w i th no n-ideal 



MHD simulations such as those of 



Fleming and Stone! (J2QDJ) and ISano and Stone! J2QQ3)- Y et 
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another possibility is that differential illumination of the disk somehow produces vorticity. Since 
the temperature of most circumstellar disks is controlled by stellar illumination, small variations in 
illumination could produce hot and cold spots in the disk that interact to produce vortices. The final 

likel y 



198' 



possibility that we consider is the generation of vorticity via barocli nic instability which i s 
to operate in disks whose vertical stratification is close to adiabatic ((Knobloch and Spruitl 
The nonlinear outcome of this instability in planetary atmospheres is the formation of vortices, 
although it is far from clear that the same outcome will occur in disks. Finally, we note that a 
residual amount of vorticity can be generated from finite-amplitude compressive perturbations. We 
have performed a series of runs with zero initial vorticity and perturbation wavelengths on the 
order of the scale height, and the results are qualitatively similar to Figure I57T1 with the shear stress 
reduced by nearly two orders of magnitude. 
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Summary and Outlook 



Angular momentum transport is key to the evolution of accretion disks. In ionized disks, momentum 
transport is likely to be mediated internally by MHD turbulence generated by the MRI. Despite 
the success of this local shear instability in elucidating the accretion process in ionized disks, the 
complexity of its nonlinear outcome has raised a whole new set of questions regarding its effects upon 
disk evolution. The answers to most of these questions will require the use of three-dimensional 
numerical simulations. I discuss some of the remaining open questions in Section [6.11 and propose 
some simple numerical experiments that will attempt to answer them. 

The mechanism driving accretion in weakly-ionized disks remains unclear. I summarize the 
main results of this dissertation in Section I6.2( results which mostly argue against a turbulent 
transport of angular momentum in these disks. I also discuss some possible directions to take in 
further pursuit of such a mechanism. 

The fact that decades of research have not uncovered a robust turbulent transport mechanism 
for weakly-ionized disks raises the possibility that at least some of these disks (or portions of them) 
are stable and do not accrete in a steady state as the standard disk model assumes. Proposals for 
exploring the implications of this possibility are discussed in Section 16.31 

Finally, while modeling turbulent shear stresses as an alpha viscosity has turned out to be 
useful phenomenologically, representing disk turbulence as an alpha viscosity has its limitations 
(see Section I6.1.2|) , and any model results that depend upon an accurate representation of this 
aspect of disk physics are therefore suspect. The fundamental understanding of turbulent shear 
stresses in disks that has begun to emerge in recent years has opened up an exciting opportunity 
for developing physically-motivated disk models based upon a first-principles treatment of disk 
turbulence. I conclude in Section IQ1 with a discussion of proposed research along these lines. 

6.1 Ionized Disks 

The phenomenological approach to modeling turbulent transport in accretion disks has been chal- 
lenged by the recent improved physical understanding of that transport in ionized disks. At the 
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same time, there are important issues with regard to MHD turbulence and its ramifications for 
disk physics that must be understood before the standard disk model can be replaced with models 
based upon a more accurate representation of turbulent transport. This section discusses a few of 
these issues and proposes ways in which they can be investigated. 

6.1.1 Transition from 2D to 3D MHD Turbulence 

An important assumption underlying standard disk theory is that the global and local disk scales 
are well separated; small-scale turbulence is assumed to average to a smoothly-varying flow on 
large scales. The validity of this assumption can be tested by measuring the power spectrum of 
MHD turbulence in a series of three-dimensional shearing-box simulations with horizontal scales 
much larger than the vertical scale, looking for a transition between two-dimensional and three- 
dimensional behavior. If there is a transition at some characteristic scale, presumably on the order 
of the disk scale height, this will confirm the standard picture as well as provide a guide for the 
scales at which the assumptions of standard disk theory can be appropriately applied in global 
models. If there is no transition to smoothly-varying two-dimensional flow, much of standard disk 
theory is invalid. 

6.1.2 Dynamics of MRI Turbulent Stresses 

Turbulence is an extremely complex phenomenon, and the standard approach of modeling it as 
a viscosity is clearly oversimplified. For example, there are key dynamica l properties of M HD 
turbu lence, such as the elastic properties that produce magnetic tension logilvie and Proctor . 



2003), that cannot be modeled with a viscosity (M cKinnev and Gammie 



2002). In addition, the 



standard disk model assumes that the turbulent shear stress is isotropic, whereas MHD turbulence 
is inherently anisotropic. Any model results that depend u pon an accurate representation of the 



turbulent shear stress are therefore suspect. 



Qgilviel (J2003J) has developed an analytic model for 



MHD turbulent stresses consisting of evolution equations for the turbulent stress tensor. Such a 
model could provide the basis for more realistic analytic and numerical studies of accretion disks, 
as well as be useful for global disk simulations since it would allow one to represent the underlying 
turbulence accurately over long time scales without having to resolve the small-scale structure. I 
will attempt to test this model against local numerical simulations of MHD turbulence and seek to 
constrain the values of its free parameters. 
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6.1.3 Interaction of Waves with Turbulence 



Another key assumption underlying many disk studies is that waves will propagate through turbu- 
lence relatively unhindered. For example, models that have been developed to explain quasi-periodic 
oscillations (QPOs) in bin ary systems often invo ke characteristic oscillations of the accretion disk 



as the source of the QPOs ([Wagoner et al 



2001). The implicit assumption in these models is that 



the disk oscillations are negligibly affected by the disk turbulence. As another example, an alytic 



1983; 



Papaloizou and Lin 



1995) have 



studies of warped accretion disks (|Papaloizou and Pringle . 
shown that warps will propagate as non-dispersive waves if the turbulence in the disk is sufficiently 
small. The turbulence in these studies is modeled as a shear viscosity, but it is not clear how the 
turbulence will interact with warp propagation self-consistently. As a final example, the standard 
disk model assumes that thermal energy is both generated and radi ated locally, which resu lts in 



a temperature dependence with radius that is not always observed ((Robinson et al 



1999). An 



alternative to this local dissipation process is a global process whereby waves are excited at one 
point in th e disk and carry en ergy and angular momentum to another point in the disk before 



dissipating ([Adams et al 



1988). Again, wave propagation through the disk is assumed to proceed 



unhindered by the turbulence. 

The general problem of the interaction of waves with turbulence has been well-studied ana- 



lytic ally, particularly in the con text of solar oscillations ((Goldreich and Kumar 



also 



Farmer anc 



Torkelsson et al 



Goldreich 



10 



1988, 



1990; see 



2004), but numerical studies along these lines are less advanced (see 



2000| for one example). I will seek to further our understanding in this area by 



investigating, via numerical experiments, the effect of turbulence on waves propagating through a 
disk. A wave or superposition of wave modes calculated from linear theory will be inserted into a 
magnetized turbulent numerical model, and the subsequent evolution will be compared with the 
linear theory results. The outcome of these experiments will determine whether or not the neglect 
of wave-turbulence interactions is valid. If there is significant interaction, these experiments may 
provide a means for developing a predictive theory of wave propagation in turbulent disks. 



6.2 Weakly-Ionized Disks 

The prospects of discovering a turbulent transport mechanism in weakly-ionized accretion disks 
appear to be dim. Since proving stability is much more difficult than demonstrating instability, 
however, the search for such a transport mechanism continues. The most promising mechanism 
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based upon the work discussed in this dissertation is the driving of compressive motions by vortices, 
but even this faces serious difficulties due to the decaying angular momentum flux and the potential 
for the vortices to be unstable in three dimensions. It would be useful to redo in three dimensions the 
simulations discussed in Chapter[51 to determine the rate at which the vortices become unstable and 
the nonlinear outcome of such an instability. Testing some of the vorticity-generation mechanisms 
discussed in fc )5.5l is also a necessary next step in understanding the relevance of this overall picture 
for driving accretion in weakly-ionized disks. 

The possibility of a bypass transition to turbulence due to the transient amplification of linear 
disturbances followed by a nonlinear feedback from trailing shwaves (shearing waves) into leading 
shwaves has not been fully explored. The effects of aliasing shown in Figure E.15l clearlv demonstrate 
that such a feedback mechanism can result in the overall growth of linear disturbances into the 
nonlinear regime. While the feedback in Figure 14.151 is entirely numerical, the high-resolution 
requirements for tracking these shwaves (as can be seen in Figure FO]) implies that the nonlinear 
outcome of transient amplification has not been tested at the resolutions we employ for the runs 
dis cussed in Chapter IH 



Vanneste et al 



(1998) calculated thrcc-shwavc interactions in the unstratified (incompressible) 
shearing sheet and found that there is feedback from trailing shwaves into leading shwaves for a small 
subset of initial shwave vectors. It would be of interest to revisit this calculation in the stratified 
shearing sheet. One key difference between the stratified and unstratified shearing-sheet models is 
that in the latter case all linear perturbations decay after their transient growth, whereas in the 
former case the density perturbation does not decay. This implies that quasilinear interactions are 
more likely to take place. A comparison of Figures 15.91 and 14.151 indicates that feedback due to 
aliasing in the unstratified sheet does not result in any overall growth, whereas feedback in the 
stratified sheet does. 



6.3 Layered Disks 



A more fruitful line of research may be to simply assume that a weakly-ionized flow is stable. 
Accretion in that case could proceed in surface layers that are ionized by non-t hermal radiativ e 



processes, with the mid-plane of the disk remaining inactive (see Figure 16.11 and 



Gammie 



1996). 



One of the few numerical studies of layered accretion has shown tha t there may be some way e 



transport from the active, MHD-turbulent zone to the inactive zone ((Fleming and Stone 



2003) 
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This work could be expanded upon with a more realistic treatment of the vertical structure and 
vertical boundary conditions, as well as a closer inspection of the stresses in the inactive zone, 
including an investigation of their dependence on simulation parameters. Understanding these 
processes would be useful for developing more sophisticated models of layered disks. 



Figure 6.1 Ionization structure of YSO disks. The inner disk is coupled to the field via thermal ion- 
ization. At larger radii the surface layers are coupled but the midplane (hashed) is inactive. At still 
larger radii the density is lower and nonthermal ionization provides effective coupling throughout 
the disk. 

In addition, since the dynamics of the dust layer in protoplanetary disks (disks which are likely 
to be weakly-ionized) have important implications for understanding planet formation, a useful and 
natural extension to the research proposed here is to incorporate gas-dust interactions in some of my 
calculations. The effects of turbulence and wave transport on the dust layer could be investigated 
by including dust particles of various sizes in numerical simulations similar to those described in 
Section 16.1.31 In addition, interactions could be calculated self-consistently by incorporating dust 
dynamics into layered disk models using a two-fluid approach. 




thermal 
ionization 




ionization 



critical radius 
- 0,1 AU 
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6.4 Advanced Physical Disk Models 



The majority o:' 


analytic ( 


Shakura and Sunvaev 




1973 




Lvnden-Bell and Prinele 




1974: F 


rinele 


1981; 


Hartmann 


1998 


) and mimerical ( Ieumenshchev and . 


4.bramowiczl. 


199?J; 


Stone et al. 




Ieumenshchev and Abramowicz. 


2000; 


Ieumenshchev et al. . 


2000; 


McKinnev and Gammie 


2002) 



accretion-disk studies incorporate the assumption of an alpha viscosity to model disk turbulence. 
As discussed in fc !6.1.2( such an assumption has its shortcomings. A fundamental understanding of 
MHD turbulent stresses can be used to develop advanced disk models that will enhance our ability 
to explain observations of accretion systems. 

Numerical models that rely on an alpha viscosity could be improved upon by incorporating a 
more sophisticated model for the turbulent stresses based upon simulations and theoretical studies 
of turbulence, such as those described in Section 16.1.21 In addition to being useful for modeling 
accretion systems, these models would provide a more reliable bridge between local and global 
MHD simulations and possibly be a means for understanding the complex results of global MHD 
simulations at a more fundamental level. 

Analytic models that incorporate more sophisticated turbulent-stress modeling could also be 
developed in an effort to improve upon the standard disk model. In addition, there may be key 
analytic results that need to be revisited with an enhanced understandi ng of turbulent stresses , 
such as the work on wave propag ation in warped disks discussed in ^16,1.31 ( Pa paloizou and Pringk , 



1983 



Papaloizou and Lin . 



1995). All of these proposals represent potential progress towards a 



first-principles understanding and modeling of accretion-disk systems. 
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The Boussinesq Approximation 



We demonstrate here that the Boussinesq approximation to the linear perturbation equations is 
formally equivalent to a short- wavelength, low-frequency limit of the full set of linear equations. 
We perform the demonstration for the stratified shearing-sheet model since the standard shearing 
sheet is recovered in the limit of zero stratification. 

Combining equations (|3.45|) through (|3.49|) into a single equation for 5v x yields the following 
differential equation, fourth-order in time: 



F 4 5v x 4) + F 3 6v® + F 2 5v x 2) + FxSvM + F Sv x = 0, 



(A.l) 



where 



y y z> 1 k x L P J y km 2 



F 3 = -2qnk x k y (k z y + hi) 1 + 



k x L P 



(A.2) 
(A.3) 



c s^x 



kl 1 + , 



k x Lp 



+ 



NJ + k 2 
klH 2 



{kl + kl) 1 



(k 2 y + k 2 z ){(k 2 + k 2 z )[l 



k x L P 



k x L P 



+ t ,2 W +l)(j + 2) 



+ k 



k 2 H 2 

2 2(q + l)(3g + 2) 
k 2 H 2 



(A.4) 



Fi = AqQc 2 k y k x 



{kl + k 2 z ) 



1 + - 



k x L P , 



1 + 



z(3g 



+ 



K 2 



iV 2 + £ 2 



2qk x L P 
a (?+!)(? 



4g/c 2 L 2 



3A: 



/fc 2 # 2 




3qk x Lp 



(A.5) 



F = c 2 fc 2 [fc 2 (iV 2 + 2g 2 ft 2 ) + fc*(J\£ + k 2 ) 



(A£ + A£) 1 



k x Lp 



+ k 



2 2(g + l)(3g + 2) 



k 2 H 2 



(A.6) 

The above expressions have been written to make the short-wavelength limit more apparent: 
all but the leading-order terms in brackets are proportional to factors of (k x Lp)~ l or (k x H) . 
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Notice also that since one expects H/Lp <C 1 for Keplerian disks with modest radial gradients, 



U « 7r4-> (A-7) 



k x L P k y Hf L P kyHr 

(for f ^ 0) and therefore the short-wavelength limit is sufficient. One needs to be careful in taking 
this limit, however, since k x = k y f goes through zero as a shwave goes from leading to trailing. 
The approximation is rigorously valid only for f ^ 0, but we have numerically integrated the full 
set of linear equations (equations (|3.45j) through Q3.49JI with k z = 0) and found good agreement 
with the Boussinesq solutions described in §4.3 for all f at sufficiently short wavelengths. 1 
With these assumptions in mind, to leading order in (Hky)^ 1 equation ([A.ljl becomes 

k x 5v x 4) - 2qQk y 5v x 3) + c 2 k x k 2 5v x 2) + AqVLc 2 s k 2 x ky 5v x ] + 

c 2 k x [k 2 (N 2 + 2q 2 fl 2 ) + k 2 (N 2 + k 2 )] 5v x = 0. (A.8) 

If we assume dt <C c s k y , the two highest-order time derivatives are of lower order and can be 
neglected (thereby eliminating the compressive shwaves) and we have 

k 2 5v x + AqQk x k y 5v x [k^N 2 + 2q 2 U 2 ) + k 2 {N 2 + k 2 )] 5v x = 0. (A.9) 

This is equivalent to equation (|3.56|) . 

Notice also that the assumption dt ~ 0(c s k y ) applied to equation (|A.8j) yields 

8v^ + c 2 k 2 6v® = (A. 10) 

to leading order in (Hky)^ 1 . This equation is of the same form as the short-wavelength limit of 
equation ()3.16j) for the compressive shwaves in the unstratified shearing sheet, confirming our claim 
that short-wavelength compressive shwaves are unchanged at leading order by stratification. 



One must start with a set of initial conditions consistent with equations 13.461 , (13.471 , 13. 5H and 13. 5211 in order 
to accurately track the incompressive-shwave solutions. In addition, suppression of the high-frequency compressive- 
shwave solutions near f = requires k y Lp > 200, which for H/Lp — 0.1 implies Hk y > 20. 
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